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ABSTRACT 


The machine processing of spatially variant multitemporal data 
such as imagery obtained at different times requires that these data 
be in geometrical registration such that the analysis processor may 
obtain the datum for a specified ground resolution element in each of 
the sets of imagery being utilized for analysis. 

Misregistration between corresponding subsets of imagery contains 
both a displacement and a geometrical distortion component, and the 
affine transformation is postulated to characterize this misregistra- 
tion between data subsets. Search techniques utilizing the moduli 
of the Fourier Transforms of these data are developed for estimating 
the coefficients of geometrical distortion components of this model. 

Following the correction of these distortion components, the dis- 
placement is located by the crosscorrelation of a template obtained 
from one set of data, termed the reference, with the second, or back- 
ground data. This template, derived for the optimum discrimination 
of the reference data embedded in the background, is determined by 
the solution of a system of equations involving the reference data 
and the covariance matrix of these data. 



The derivation of the optimum filter includes constraints such that 
the maximum filter output, corresponding to the correct superposition 
of the reference template on the background data, is unity and the energy 
in the filter is finite. The filter obtained in this development is 
linear although it may involve a parameter requiring the solution of a 
nonlinear equation. 

The performance of the crosscorrelation algorithm is evaluated 
using ideal data obtained by convolving an array of computer generated 
random numbers with a two-dimensional lowpass filter having a specified 
impulse response. The results obtained from these data generally 
substantiate the conclusions drawn from the analysis of this algorithm. 
The correlator output is then obtained for noise free and distortionless 
line scanner data. In these data the reference is selected as a subimage 
of the background data, and the data are selected to typify line scanner 
imagery. Multitemporal data are processed with the algorithms developed 
for the noise -free data to evaluate the applicability of this filter 
to the conjugate point problem. 

It is demonstrated that the crosscorrelation of the template 
derived from the reference data will not yield useful results unless 
the geometrical correction of the data is implemented. The Fourier 
transform search techniques are used to estimate the distortion model 
coefficients, and a bilinear interpolation algorithm is utilized to 
correct the imagery. Results of the processor output using the 
corrected data are given. It is shown that the optimum filter yields 
a more discriminable peak of the correlation surface at the correct 
superposition of the reference template on the background than does the 



filter chosen as a sub image of the reference data Itself. 
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CHAPTER 1 
INTRODUCTION 

1.0 Processing of Remotely Sensed Data 

The large scale application of remote sensing techniques to the 
monitoring and evaluation of the environment is rapidly becoming a 
reality. This realization is due to many technical and scientific 
advances resulting in increasing availability of multlspectral photo- 
metric and radiometric instruments, aircraft and satellite instrument 
platforms and large data processing systems [28,42], The data 
obtained by these multlspectral instruments are images representing 
the spatial, spectral, and temporal characteristics of the area 
under investigation in contiguous wavelength intervals throughout 
the visible region of the electromagnetic spectrum as well as 
selected wavelength bands in the infrared and microwave portions 
of this spectrum. 

A block diagram of a multlspectral remote sensing instrumentation 
system is illustrated in Fig. 1-1. To analyze the large quantity of 
data made available by such a source, statistical pattern recognition 
algorithms [28] have been developed, and such pattern recognition 
procedures are illustrated in Fig. 1-2. The feature extractor selects 
from the pattern those attributes which "separate" a particular 
pattern as well as possible from the set of all other patterns. The 
output of the extractor is a k-tuple and is usually of lower dimension 






3 








4 

than the pattern vector. The pattern classifier then assigns the 
input vector to one of r-classes based on some predetermined decision 
rule. 

An implicit assumption in such processors is that any one input 
vector is derived from a common ground resolution element. This 
requirement is usually satisfied if all the elements of the n-tuple 
are taken at time t ■» t^ with an instrument employing a single optical 
path and negligible differential optical and electrical delays are 
introduced into the channels following the dispersion of the input 
optical signal. 

There has been recent interest in extending the pattern classifier 
input vector to include a multitemporal variation. However, such data 
sets are generally misregistered because after obtaining the first data 
set, the instrument platform cannot be made to follow the same path 
to within a resolution element of the instrumentation system. The 
misregistration is a dynamic quantity and can be considered quasi- 
static only over that length of data defined by the dynamics of the 
instrument platform and the conditions of the atmosphere in the case 
of aircraft platforms. Therefore, for multitemporal data, a pre- 
processing operation must be implemented for removing this misregistra- 
tion before the data analysis system can address the same spatial 
data element in each of the data sets being used for analysis. This 
requirement for addressing a data element in the various data sets 
is illustrated in Fig, 1-3. 

1.1 Multi temporal Data Registration 

A multi temporal data registration system, shown in Fig. 1-4, can 




Overlay Processor 
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be divided into two relatively distinct parts with this dichotomy 
being a functional division. The ultimate objective of an overlay 
system, given two sets of imagery A and B, is to process set B such 
that its image under the transformation T is in geometrical registra- 
tion with set A. Thus the first part of the overlay problem is given 
enough Information to define the transformation T, how can the mapping 
be efficiently implemented for the rectification of set B under the 
constraints of computer memory size and data through-put requirements? 
This correction of a misregistered data set is illustrated in Fig. 1-5 

The solution of this part of the problem usually Involves the 
division of the total data set into smaller Images, or blocks, which 
are more readily handled. A transformation is derived for each of 
these blocks and the rectification is carried out. If it is required 
to have continuity of the image across the boundaries of the blocks, 
this condition can be approximated by making each block sufficiently 
small such that the differential scale factor or rotation will not 
alter the block size more than some fractional part of the size of 
an image resolution element. The alternative procedure of formulat- 
ing the problem such that the blocks of set B under the transformation 
T are constrained to be continuous across the boundaries may be 
possible. 

The second part of the data overlay problem is the task of 
determining the transformation T relating any two blocks of data. 
Assuming that a model has been determined for T, the parameters 
of this model must be estimated from the data since the misregistra- 
tion is a pairwise property of the data sets and is implicitly defined 



Distorted Image Rectified Image 


Figure 1-5 Rectification of Distorted Imagery 
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in these data. 

For the unconstrained problem, a conanonly used model for 
representing the misregistration between data blocks is the two- 
dimensional polynomial 

H-i N-i . 

2 - E 2 a,, xj x; i-1,2,. (1-1) 

j-0 k*0 J 1 

The principal advantage of this polynomial model is that equation (1-1) 
is linear in the coefficients a^. If a set of corresponding, or 
"conjugate" , points is known for each of the blocks constituting a 
set of Imagery, the coefficients of the transformation for each of 
these blocks are readily obtained by use of a least- squares procedure. 
Thus the requirement arises for determining conjugate points in sets 
of multitemporal imagery. 

Signal processing techniques used to determine conjugate points 
almost always use correlation between data sub-blocks as a measure 
of the similarity of these data subsets. To illustrate this 
processing, assume the data to be a discrete set obtained by sampling 
a continuous scene on a rectangular grid and the reference data 
subset S c A to be selected from data set A. A background data sub- 

— T 

block is chosen from data set B and the following inclusion relation 
is assumed 

T 

s, -> s; c*b (1 ‘ 2) 

where T is the transformation relating the data blocks. 

The correlation surface 

C(u,v) - E l « r (i,j) \(u + i, v+ j), s r 6 ® b € S b (1-3) 



1 ° 

is confuted for the set of indices (i, j), and the peak of this 
surface is assumed to be the location of the correct superposition 
of on £5^* The conjugate point pair is then taken to be the 
location of j> r and the corresponding subset of in their respective 
coordinate systems for this maximum of the correlation surface. 

1.2 Objective of this Investigation 

This investigation is addressed to the problem of conjugate 
point determination in multitemporal imagery in which the task is 
viewed as a problem in signal theory. With conjugate points defined 
as in section 1.1, this problem can then be phrased as "given a 
reference set from data set A, what signal processing operations 
will optimally discriminate the corresponding subimage in data 
set B?" This problem is treated in the sequel with consideration 
being given to both the noise and misregistration which exist 
between the data sets. 

1.3 Previous Investigations 

The implementation of an image registration system on a general 
purpose computer is a relatively new field although the roots of 
this task extend into antiquity. Only recently has there appeared 
in the literature any significant amount of work on the image 
registration problem. The work which has appeared as well as the 
experimental effort in this study have utilized digital computers. 

In measuring the correlation surface, Anuta [4] used the 
normalized correlation coefficient, computing quanlty 
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M N 

£ £ s (i + u , j + v) s (i,j) 

R(u,v) - 1 ~ 1 (1-4) 

r M N 1/2 

£ £ s‘(i + u, j+v) • SQ 

L i-1 j-1 b J 


where 


SQ - £ £ s‘(i,j) 

1-1 j-1 


with s r € and € S^. In a later paper Anuta [5] utilized the 
Fast Fourier Transform algorithm for increasing the efficiency in 
evaluating equation (1-4). 

Rather than use the normalized correlation coefficient as a 
measure of image similarity, Baraea and Silverman [7] used an 
absolute difference metric for their similarity algorithm 


e - £ £ 

i J 


l* r (i»j) - s b (i + «» J + v) 1 


(1-5) 


The rational for this algorithm is that if the reference data S^. 
are not near the correct superposition on S^, the error will Increase 


rapidly and a threshold can be established such that if this boundary 
is exceeded, the checking for similarity for this observation is 
terminated. It is suggested that the average number of operations 
for testing for the correct superposition can be significantly 
reduced compared to the correlation algorithm. 

Arcese et al . [6] developed a filter which maximized the ratio 
of the square of the filter output at the correct juxtaposition of 
S on S.. Their filter was developed for the noise free case and 



12 


with no constraints on the filter characteristics. In addition 

the data were assumed to be from a process having a simple exponential 

covariance function. 

The most complete system for the overlay of large data sets was 
described by Lillestrand [46]. This work incorporates some of the 
ideas that were independently discovered during this investigation; 
particularly in the use of pre-processing to estimate and reduce the 
geometric distortion so as to increase the similarity of corresponding 
data sets. Hall et al. [33] discuss a method of registering cloud 
photographs. Their procedure employs a crosscorrelation method 
similar to that of Anuta for locating corresponding geographical 
features in the two data sets. 

1.4 Outline of Investigation 

The misregistration of the imagery is discussed in Chapter 2. 

It is observed that the geometrical components of the misregistration 
between these data can be considered to have been introduced by a 
spatially variant operator. In view of the analytical difficulties 
in formulating a general solution to this problem, a sequential method 
is proposed in which the geometrical distortion components are 
estimated and removed prior to determining the displacement components. 
It is postulated that the affine model adequately represents, for the 
purposes of this investigations, the misregistration between data 
subblocks; an algorithm utilizing the modulus of the two-dimensional 
Fourier transform of these data is developed for estimating the 
geometrical components of this model. 
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The optimum filter for discriminating the data set embedded 

in the background data is treated in Chapter 3. This filter is 

derived under constraints on both the filter output and the filter 

energy. Results obtained in this chapter include expressions for 

the filter for both noise free inputs and for inputs assumed to 

contain additive noise with statistical characteristics given by 

the covariance matrix K . 

— n 

In Chapter 4 experimental results are presented using both ideal 
data obtained by convolving computer generated random numbers with 
lowpass filters having known Impulse responses and Imagery obtained 
with a line scanning instrumentation system. Plots of the correlator 
outputs for various input data sets of interest are given, and the 
processing gains of the correlator are given for both ideal and 
multi temporal line scanner data. 

A summary of the results obtained in this investigation is 
presented in Chapter 5. 



CHAPTER 2 


DATA MISREGISTRATION 

2.0 Introduction 

Misregistration of data taken at different times occurs when, 
after having obtained the first set of data, it is not possible to 
control the positioning of the instrument platform such that it will 
repeat the previous path to within the dimensions of the spatial 
resolution of the instrumentation system. For certain kinds of 
processing of multitemporal imagery, it is necessary to Implement 
a data pre-processing operation such that the data are placed into 
registration. Following such an operation on the data, the computer 
or other processor is able to address (or obtain) the data from a 
common ground resolution element in the several data channels being 
used for analysis. 

It is the purpose of this chapter to develop an appropriate model 
of the misregistration process and to investigate some methods of 
evaluating parameters of such a model. The first topics discussed 
are the characterization of the data, and definitions and nomenclature 
that are used subsequently. Models for regional misregistration are 
then discussed, and the affine transformation is proposed as a suit- 
able representation of the misregistration process. Processors for 
estimating the misregistration are introduced and it is concluded 
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that for the data of interest in this study, the geometrical 
distortion and displacement parameters can be determined sequentially. 
The chapter concludes with the development of Fourier transform 
methods for estimating the coefficients of the affine model which 
represent the geometrical distortion components of the misregistration 
between an image pair. 


2.1 Data Source 

The data source s(x) is assumed to be a real, stationary, wave- 
number limited, two-dimensional random field. Each realization of this 
process, is shown in Fig. 1-1, is obtained by impulse sampling of s(x) 
on a bounded rectangular region R. Thus the sampled image is 


- “ Wp “ •(£)•' * "(*) 


( 2 - 1 ) 


where tt Is the two-dimensional finite sequence of Impulses spaced at 
unit distance 


tt(x) » S 2 6(x -n , x -n„) rect 

n^»-°° n2“-°° 1 2 


(¥*) (¥*) 


( 2 - 2 ) 


with rect <£) - {J; 


and (P 0 >4 0 ) c ^e center of the data aperture. 


The set 


S - {s ij : 1 < i < p, 1 < j < q) 


may be represented by 


the column vector 




(2-3a) 


s 



where the elements s, are the column vectors 
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2 * 


(2-3b) 


iq J 


These data are also quantized; for each sample s^ € S, the value Is 
a D a^, 0 < k < K-l where the a^ are the quantization levels and K 
Is the total number of levels provided by the quantizer. 

Each sampled image S> obtained from s(x) is called a quantized 
picture function, or synonymously a digital picture function. Each 
picture function S is an element in the pq-dimensional Euclidean 
space, and the distance between any two digital picture functions P 
and £ is then 


d(p,q) - ||P all’ - (P.,2) 


1/2 


(2-4) 


[(p - 2 ) t (p - 2)] 


1/2 


{ [(£■ - 2 > <1 * 9 ) T ] } 


1/2 


where tr (*) is the trace of the matrix. It should be noted that once 
a digital picture is obtained, it can be regarded as a known 
deterministic, discrete, two-dimensional signal. However, the 
statistical properties of this picture, such as variance or bandwidth, 
are in fact related to these same properties of s(x). 

Considering now the content of the picture function, an "object" 
is a subset of S which represents some identifiable physical entity in 
the object plane. The distance between two objects in a picture 
function is defined as the sum of their row and column distances. 
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Thus for two objects located at (x^, x 2 ) and (y^, y 2 ) respectively 

A row " l x 2 ' y l\ < 2 ‘ 5a > 

L col " K • y l\ (2_5b) 

and 

<«*. 2> - 4 ro „ + A eol (2- 6 > 

A "scene" is composed of a collection of objects together with 
their mutual geometrical relationships; the term scene Is synonymous 
with the term Imagery If the data are assumed continuous, whereas It 
Is synonymous with digital picture function when the data are in a 
discrete format. 

Geometrical distortion between two scenes r(x) and s(x) Is defined 
to be the mapping D of one scene onto the other. Such a mapping can be 
expressed as 

r(x) - D 0<x)] - s[5(x), Tl(x)] (2-7) 

where £ , T| are real valued functions and this mapping Is one-to-one 
and onto. Distortion manifests itself In the differences In the 
distances between corresponding objects in the two digital pictures 
R, S>, where 

R = r(x) • 2 tt(x) (2- 8a) 

S - s(x) • 2 tt(x) (2- 8b) 

Translation is defined as 

r(x) - T[s(x)] - s(x^ + tj, x 2 + t 2 ) (2-9) 
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Misregistration in general can be expressed by the product of the 
operators D and T as 

r(x) - (TD) s(x) - T{D[s(x)]} (2-10) 

which for digital picture functions results in 

r(x) • 2 tt(x) - (TD[s(x)]} • 2 tt(x) (2-11) 

In the sequel the discussion often requires a qualification as to 
the size of the data set being considered. To this end the imprecisely 
defined terms local, global and regional will be used. The term local 
is used to indicate that the data set consists of an element s^ 
together with those immediately surrounding elements (s^) such that 

d (8 tj , 8 k i) - |i - fc | + | j - X| < N (2-12) 

with N as integer. 

The term global is used In the sense that the entire data set is 
being considered. Between these extremes, the term regional is used 
to define a data set for which one choice of distortion model 
coefficients serves to adequately represent the misregistration between 
the data and its conjugate image. The size of a region is very dependent 
on both the assumed misregistration model and on the data itself. 

2.2 A Regional Misregistration Model 

Remotely sensed images typical of those obtained by a line 
scanning instrument are shown in Fig. 2-1. As illustrated by this 
figure, seldom is the misregistration between two data sets 
characterized solely by the displacement of one image with respect 





Example A 


Example B 


Figure 2-1 Multitemporal Imagery Exhibiting Misregistration 
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to the other. Rather, a geometrical component, as defined by 
equation (2-7) ia also present. A commonly used model for character- 
izing this misregistration is the two-dimensional polynomial 


y 


i 


N-l N-l 
E E 
j-0 k-0 



k 

X 2 


1 * 1,2 


(2-13) 


where x and are respectively the coordinate systems of the reference 
and background data sets. 

The utility of the polynomial model lies in the linearity of 
equation (2-13) as a function of the coefficients a^. Thus, this 
model is useful in least-squares procedures where a set of correspond- 
ing image points is given and an equation of a given order which best 
relates these points is desired. Table 2-1 lists several polynomials 
obtained from equation (2-13) for degree N-l less than or equal to 
three. Also given is the number of coefficients *8^, i - 1,2 which 
must be specified for a polynomial of a given degree. It is evident 
that the number of coefficients required for a model increases rapidly 
with the degree of the polynomial. 

As a function of the coordinate variables x, equation (2-13) is 

linear only if the coefficients *a^ t 0, with i ■ 1,2. The terms ^a^ 

2 

and a„ represent the displacement of the origins of the two 
coordinate systems x and This nonlinearity is removable by the 
coordinate translation 


y i ’ x i " a 00 ’ 1 “ 1 » 2 


(2-14) 


However, the a^ are unknown and are the quantities of ultimate 


interest. For all indices j ,k >_ 2, equation (2-13) describes a 
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Table 2-1. No. of Parameters Required for Two- Dimensional 
Polynomial Misregistration Models. 


Common 

Name 

Degree 

N-l 

No. Parameters 
Required to 
Determine Model 

displacement 

0 

2 

linear 

1 

4 

affine 

1 

6 

projection 

2 

8 

quadratic 

2 

12 

cubic 

3 

20 
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nonlinear transformation. 

The data of interest in this study are obtained using a 
mechanical line scanning instrument mounted in an aircraft and 
operated at an altitude of approximately 5000 feet. As illustrated 
in Fig. 2.2, this imagery is obtained by recording the scene radiance 
from an effective ground resolution element as this resolution element 
is moved along an approximately linear locus by the rotation of a 
mirror within the instrument. Adjacent scan lines are spaced by the 
motion of the aircraft, and the angle a between the direction of 
the scan lines and the aircraft ground track is 

a - Z + y ( 2_15) 

where y is the crab angle of the aircraft with respect to the ground 
track. 

In this study attention is directed toward the misregistration 
of sub images of data, termed regions, which are used for determining 
geometrically corresponding points in the two data sets. The size 
of a region ranges between 32 x 32 and 128 x 128 picture elements. 

With data sets of this size and with any changes in the aircraft 
attitude necessarily limited by its dynamics, it is assumed that 
each region of data consists of a sequence of equally spaced, linear 
scans of the scene. Each scan line is at an angle 

» = ? + Y„ (2-16) 

o z o 

with respect to the aircraft ground track and the ground track has the 
angle p with respect to a pre-assigned coordinate system. These 
quantities are illustrated in Fig. 2-3. 
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Figure 2-2 Line Scanning Geometry 
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It 18 observed that the regional misregistration then has the 
following four components: 

(1) scale 

(2) rotation 

(3) skew 

(4) displacement. 

These misregistration components are assumed to be characterized 
by the affine transformation 


_y»Ax + _t«Bx 


where 


(2-17) 



V 


■ V 


V 

z - 


, * " 


. 1 “ 



■ y 2. 


- x 2 


- t 2- 


and 


A 


*11 *12 
■ *21 * 22 - 


is a nonsingular matrix characterizing the geometrical components of 
the misregistration. The quantities are illustrated in Fig. 2-3, 
where x and y_ are the coordinate systems and 


t - 


X 10 - *10 
X 20 " y 20 


(2-18) 


is the displacement of some known point in each data set. 
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It Is illustrative to examine the distortion matrix A in detail, 
identifying those coefficients which reflect the various distortions. 
For scaling differences along directions parallel to the coordinate 
axes, the distortion matrix is 


A - 



0 


0 


22 


(2-19) 


where the a^ are the scaling coefficients. For rotation, A becomes 
the orthogonal matrix 


where 


COS P 

sin P 

-sin P 

C 08 P 

-P 2 



(2-20) 


is the difference in angular orientation of the ground tracks of 
each data set with respect to its coordinate system. 

Skew distortion is introduced by a distortion matrix of the form 


A - 


1 1 



( 2 - 21 ) 


and it is readily verified that the angle a, defined by equation (2-16a) 
is 


” «i - <*2 " tan _1 (a 21 ) 


a 


(2-22) 
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The distortion matrix A is some combination of these components 
distortions. The necessity to determine the elements of A from the 
two regional data sets and a procedure for accomplishing this task 
for certain classes of data are discussed in the subsequent sections. 

The significance of the affine model assumption is that the mis- 
registration can be interpreted as consisting of two components; the 
first is the displacement of corresponding regions and the second is 
the characterization of the geometrical distortion by a linear model. 

It will be established in the sequel that for certain classes of data 
and under very reasonable assumptions, Fourier transform methods may 
be used to separate these two misregistration components and straight- 
forward search techniques are available for estimating these distortion 
parameters in the spatial frequency domain. 

2.3 Misregistration Processor 

The function of the misregistration processor is to identify 
corresponding points in the two sets of imagery being processed. 

Since these data are ordered by a Cartesian coordinate system, this 
identification requires the determination of the translational com- 
ponent t of the assumed misregistration model. 

The general problem of conjugate point identification is 
Illustrated by the block diagram of Fig. 2-4. The operator D 
introduces the misregistration; the noise N is assumed to be additive, 
independent of S and with a constant power spectral density. The 

filter H is to be determined such that with the reference data set S 
— — r 

selected from j>, the filter output when overlays its conjugate data 
set in £5^ is maximally discriminate from the outputs at all other 
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N 



Figure 2-h Block Diagram of a General Misregistration Processor 
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spatial juxtapositions of on . 

An interesting aspect of this problem is developed by examining 
the operator D. For this operator to be spatially invariant it must 
have the property of commuting with the translation operator T [56]. 
Assuming continuous operators on f(x) 

D[T f(x)] - T [D f (x) ] (2-23) 

where T is the translation operator 

T[f(x)] - f^ + t^ x 2 + t 2 ) (2-24) 

The operator D induces the mapping 

D[f(x)] - f[5(x), Tl(x)] (2-25) 

where 71 are real and the mapping is one-to-one and onto. Applying 
the criterion of equation (2-22) to the scene s(x) gives 


D{T[ s(x)]} - s[5(x - t), TKx - t)] (2-26) 

T{D[s(x)]] - s[5(x) - t x , T1(x) - t 2 ] (2-27) 

and it immediately follows from equations (2-26) and (2-27) that 

D(T[s(x)]) 4 T{d[ s(x) ]} (2-28) 

From equation (2-28) it is concluded that 

2 tt(x) • D[t[ s(x)]} 4 2 tt(x) • T[D[s(x)]] (2-29) 

and D is a spatially-varlant operator. 


In view of the fact that the transformation D is spatially variant, 
the determination of H for the general misregistration is not pursued 


30 


farther. Rather a sequential procedure is developed, illustrated 
by Fig. 2-5, where the geometrical distortion components are estimated 
and removed before the signal is processed by the filter H. 

For the class of data of interest in this study, additional 
simplifying assumptions are made for purposes of formulating the 
distortion estimation problem in a straightforward manner. These 
assumptions are as follows: 

(1) neglect the earth's curvature 

(2) neglect scan-line foreshortening 

(3) assume uniform illumination 

With the affine model assumed to describe the misregistration and these 
simplifying assumptions, Fourier transform techniques for estimating 
the geometrical distortion are discussed in the subsequent sections 
of this chapter. 

2.4 The Two-Dimensional Fourier Transform under an Affine 

Transformation . 

Under the assumption of an affine transformation as the model 
representing the distortion between two images, the expression re- 
lating the spatial frequency domains of these images is readily derived 
using the two-dimensional Fourier transform. The significance of 
this expression is that under some very reasonable assumptions, 
estimates of the distortion coefficients can be readily made for the 
class of data of interest in this study. The analysis of this section 
will use continuous functions, as no generality is lost and the 


notation is much clearer. 
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The two-dimensional Fourier transform is defined as 

CD 


F(u) 


J J £( i> 

— CO 


exp £-J2 tt(u,x)J dx 


(2-30) 


and in polar coordinates 
® 2tt 

G(p,0) - J Jr g(r ,9) exp ^-j2TTrp cos(9 - 0)J dr d9 (2-31) 

o o 

It is shown in Appendix A that if there exists a function A(x) such 
that 


4(x) - D[f(x)] - f[B x] - f(£) 


(2-32) 


where B is an affine transformation, then 

F(v) " JJ|" ex P [-j2n(u,t)J F |J (A X ) T a ] (2-33) 

where 


dy l 

^1 

Bx l 

9x 2 

^2 

*2 

9Xj 

Sx 2 


is the Jacobian of the transformation. It follows from (2-33) that the 
modulus of the transform is independent of the displacement vector t. 


M(v) - j^F(v) 



■yjy {|f[(a _1 ) t u ]| 2 } 


(2-34) 
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With proper interpretation of (2-33) and (2-34), several interesting 
properties of the two-dimensional Fourier transform can be emphasized. 

The Fourier transform is a decomposition of f(x) into a linear 
combination of basis functions of the form 

exp ^-j2rr (u,x)J - exp £-J2 tt(u 1 x x + u £ x 2 )J (2-35) 

For a particular pair of coefficients (u^,u 2 ), the corresponding 
elementary function has zero phase along the line described by 
u 

x 0 - x. + — , n is integer (2-36) 

2 u 2 1 u^ ° 

and the wavefront has the direction 
-1 U 1 

9 - tan 1 — (2-37) 

u 2 

From Fig. 2-6 it is easily seen that the spatial period of the wave- 
front is 


L ■ — cos 9 _ 

u. 


K * - 


n 


172 


(2-38) 


Thus the corresponding spatial frequency is 

, 1/2 




(2-39) 


From (2-34) the relation between coordinates systems u and v 


is 


,.-lvT 

- (A ) u 


v 


(2-40) 
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It Is illustrative to examine the effects of various geometrical distor- 
tions in the spatial and spatial-frequency domains. For scaling changes 
the distortion matrix is 

(2-41) 



where the scaling is assumed to be along the coordinate axes. The 
relationship between the spatial frequency domain coordinates is then 

-1 T 

v - (A ) u - 

This relation is illustrated in Fig. 2-7 where the dimensions of the 
rectangular blocks in the spatial domain have a width which is two 
times the height. 

For rotation the distortion matrix is 


r JL 

a ll 
0 


*99 J 


U 


(2-42) 


A - 


cos 9 


sin 9 


L-sin 9 cos 9j 

which is an orthogonal matrix. Thus 


(2-43) 


v - (A *■)■*■ u ■ A u (2-44) 

and the transform of the distorted data is also rotated by the angle 9; 
this distortion is illustrated by the modulus of the transform of the 
ideal data shown in Fig. 2-8. 



Block Pattern 


Figure 2- 



Modulus of Fourier Transform 


7 Modulus of Fourier Transform of Block Pattern 




Rotated Block Pattern 



Modulus of Fourier Transform 


Figure 2-8 Modulus of Fourier Transform of Rotated Block Pattern 






Block Pattern with Skew 


Figure 2- 



Modulus of Fourier Transform 


9 Modulus of Fourier Transform of Skewed Block Pattern 
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geometric distortion matrix becomes 

(2-45) 

domain the coordinates are related by 

u (2-46) 

This distortion and its associated modulus of the Fourier transform 
are shown in Fig. 2-9. 

The symmetry property of the modulus is the remaining property of 
the two-dimensional Fourier transform needed for this study. For real 
signals it readily follows from equation (2-30) that 

F(-u) - F*(u) (2-47) 

where the asterisk denotes the complex conjugate. Thus the modulus, 
equation (2-34) , is symmetric about the origin, and the frequency 
domain search techniques need to use only one-half of the modulus. 

2 . 5 Transform Techniques for the Correction of Regional Geometric 
Distortion . 

For data sets with the affine model characterizing the misregistra- 
tion between corresponding regions, Fourier transforms of these data 
may be used to advantage in determining the coefficients of the 
geometrical distortion components. These advantages are a result of 
both the property of the Fourier transform exhibiting the spatial 


For skew distortion the 


1 0 


L a 21 1 


and in the spatial frequency 


21 
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characteristics existing in the spatial domain and of the organization 
of these characteristics in the spatial frequency domain. The zero 
spatial frequency component is mapped into the origin of the trans- 
form domain and higher frequency components are mapped into locations 
proportional to both the value of their spatial frequency and in a 
direction from the origin characteristic of the orientation of the 
component in the spatial domain. 

With the assumption of the affine misregistration model, equation 
(2-33) gives the relation between the spatial frequency domain 
representations of the data sets. As defined by equation (2-34) the 
moduli of the two-dimensional Fourier transforms of corresponding 
regions of data are invariant under the coordinate shift _t; thus the 
coordinate systems can be chosen arbitrarily. It is further assumed 
that for data of a reasonably homogenous composition, small shifts Ad 
of the data aperture will yield moduli which can be assumed to be 
unchanged for purposes of this study. 

The differences between the moduli of the transforms of two 
corresponding regions of data provide all the information required for 
determining values for the coefficients of the linear distortion matrix. 
For the data sets employed in this study the scale factors can be esti- 
mated accurately from the flight information available. Therefore, 
only distortion due to skew and rotation must be estimated from the 
data sets. 

The computation of the Fourier transforms for the data used in 
this study is carried out using the Fast Fourier Transform technique 
[5,9] since the data are in a discrete format. 
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2.5.1. Search Algorithms for Structured Data 

The class of agricultural imagery, which is the data of principal 
Interest in this study, typically consists of a collection of 
rectangular fields with each of these fields having essentially a 
homogenous ground cover. The data usually, but not always, are 
recorded such that the flight path is closely parallel to the 
orientation of these fields. It has been observed that the modulus 
of the two-dimensional Fourier transform of data sets from this class 
of imagery typically exhibits a simple structure with the property 
that a majority of the energy in the spatial frequency domain is 
concentrated along linear loci or rays perpendicular to these field 
boundaries. The moduli of transforms typical of those obtained for 
agriculturally related imagery are shown in Fig. 2-10(a). 

The rotative and skew components of the misregistration are 
obtained by utilizing the geometrical differences between the moduli 
of the Fourier transforms of corresponding regions. This is 
accomplished by determining the angle between the corresponding loci 
of energy which are characteristic of this class of data. Assuming 
that the location of these loci can be determined, the distortion 
components are evaluated in a straightforward manner. 

(1) Rotation 

From equation (2-44) it is seen that rotation in the spatial 
domain results in the rigid body rotation of the modulus of the 
Fourier transform. Only the angle between the corresponding loci of 
the distorted and reference data transforms is required to determine 
the rotational distortion; this is shown in Figure 2-11. Thus 



(a) Structured Data Sets 




(b) Unstructured Data Sets 

Figure 2-10 Examples of the Modulus of the Two-Dimensional 

Fourier Transform 
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Figure 2-11 Location of Loci of Energy in the Modulus of the 

Two-Dimensional Fourier Transform, Rotational Distortion 



44 


the coefficients of equation (2-43) are 


a^ - &22 = cos ® 


a 12 " ~ a 21 “ — sin ® 


(2-48) 


where the sign (+) is determined by the sense of the rotation. 

(2) Skew 

The angular quantities related to the skew distortion component 
are shown in Fig. 2-12. It is assumed that a locus of energy in the 
modulus of the Fourier transform of the reference data is at an 
angle 9 with respect to the u^ axis. The corresponding locus of 
energy of the distorted data is at an angle 0 with respect to the 
reference data locus. It then follows from equation (2-26) that 
these angular quantities and the skew coefficient a^^ are related 
by 


0 - tan 


Mr* 


tan 9 


tan 9 


To determine a 2 ^, take the limit of 0 as 9 -» tt/ 2; thus 


(2-49) 


lim tan * (• 
9 -» tt/2 


1/tan 


9 + a 12 ) 


(2-50) 


- tan 


1 (— ) 
V 


The skew coefficient is obtained from the expression 


a„, - 


21 tan 0 


IT 

2 


(2-51) 
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Figure 2-12 Location of Loci of Energy in the Modulus of the 

Two-Dimensional Fourier Transform, Skew Distortion 
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(3) Rotation and Skew 

In the more general case of both a rotative and skew distortion 
component in the data, there must be at least two linear loci in each 
modulus of the transforms of the data sets to determine these 
distortion parameters. For agriculturally related data this require- 
ment is usually satisfied with the two loci being separated by 
approximately tt/2 radians. 

With reference to Fig. 2-13 assume, without loss of generality, 
that one of the loci of the reference data transform is coincident 
with the Uj axis. Thus from equation (2-49) 


1 ( 

\ 1 + 


tan 6^ 
a 21 tM "l 


) + Oj - Qj (2-52) 
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Figure 2-I3 Location of Loci of Energy in the Modulus of the 
Ttoo-Dlmentsional Fourier Transform, Rotational 
and Skew Distortion 


48 


For rectangular fields for which © 2 - rr/2 


a 


21 


1 

tan (0 2 - 0 t ) 


1 

tan 0 2 


(2-56) 


where 0 2 - © 2 - 0^ 

The distortion matrix is then the product, with attention given 
to the proper order of the matrices 


A ■ 


s 


L 21 


1 0 
1 


11 

r 12 

21 

r 22 - 


(2- 57a) 


11 


12 


S 21 r ll + r 21 ®21 r l2 + r 22 


(2-57b) 


which yields 


cos 0, 


+ sin 0. 


a 2l cos 0 ^ + sin 0 ^ + a 21 sin 0 ^ + cos 0 ^ 


(2-58) 


The algorithm Implemented in this study, for locating the loci 
containing the significant fraction of the energy in the spatial 
frequency domain, utilizes a two part procedure. The concept of the 
first operation is shown in Fig. 2-14(a). The wedge aperture is scanned 
through a total angle T in steps equal to the angular resolution of 
the aperture AY. The output of this spatial filter is the summation 
of those data falling within the region (P 2 - P^)AY. The output is 
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9 , Total Angle for 



(a) Wedge Aperture 


Maximum Value 



Figure 2-14 Use of Wedge Aperture to Locate Loci in the Modulus 
of the Fourier Transform of Structured Data 
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V i +AV p 2 


l M(p, Y^, 


i - 1, 


(2-59) 


where H(p,v) is the modulus of the two-dimensional Fourier transform 
expressed in polar coordinates and N is the number of resolvable 
positions of the aperture within the total scan angle. 

This summation necessarily involves data Interpolation since 
the data are located on a discrete grid. The w^ are nonnegative 
and the largest w^ is assumed to be obtained from the segment con- 
taining the loci being sought. 

The purpose of the second step is to locate the loci within this 


segment. A sequence of rectangular segments is selected from the 


data along the radial axis of the wedge giving the maximum output 


as shown in Fig. 2-140>). The largest data point in each of these 
segments is found, and the slope P of the linear regression line 
through the coordinates of this set of maxima is computed. Thus 


(u U • u l> <u 2i ' u 2> 


" - 2 
I (u " u i> 
i«l 


(2-60) 


U ! " M £ U li 


U 2 " M ^ U 2i 
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arid M is the number of points selected for the regression procedure. 
The estimated angle 9 between this locus and the axis Is then 

9 - tan" 1 P (2-61) 

2.5.2. Search Algorithms for Data with Complex Structure 

The class of scenes of natural or uncultured terrlan has Fourier 
transforms whose moduli have a more complex structure than the class 
considered In the previous section. The moduli of transforms typical 
of these obtained for this class of data are shown In Fig. 2-10(b). 

To identify the distortion parameters relating the transforms of two 
corresponding regions of data from this class of Imagery necessarily 
involves a search procedure. 

A possible search algorithm is to examine the modulus of the 
transform of the reference data along some locus and then search 
for the corresponding locus In the distorted data such that the 
difference between the outputs of these two spatial filters is 
minimized. A constant frequency search aperture, shown in Fig. 2-15, 
for the search locus in the transform of the reference data is 
appealing because it is both simply implemented and its corresponding 
locus under a linear transformation is readily related to the 
distortion matrix parameters. 

The output of the constant frequency spatial filter is the 

sequence of numbers w^ 

0j^ + A0 P£ 

- L 2 MCp.0^, i-'l N (2-62) 
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where M(p,0) is the modulus of the two-dimensional Fourier transform 
expressed in polar coordinates and N is the number of resolvable 
positions of the filter with the rr radian search extent. 

(1) Rotation 

For rotation the search procedure is particularly simple 
because from equation (2-44) it is seen that a rotational distortion 
introduces only a rigid body rotation of the transform modulus. 

Thus the function to be minimized is 

N , 

e - E | w ri (0) - w bi (0 + 0 o ) T (2-63) 

1 = 1 


where the subscripts r and b refer to the reference and background 
data sets repectively, and 0 q is the variable for which the minimum 
is sought. It then follows that the parameters of the distortion 
matrix A are 

a Il ’ a 22 ' coa 0 o 

a 12 “ " a 21 " - 8in (2-64) 

where the sign is dependent on the sense of rotation. 

(2) Skew 

For a constant spatial frequency search locus in the reference 
data, the corresponding locus in the distorted data domain under a 
skew transformation is an ellipse. This is shown in the following 
development. Let the circular locus be given by 


u 


2 

1 



2 

u 

o 


(2-65a) 
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From equation (2-46) the transformation of coordinates Is given by 


U 1 “ V 1 + *21 v 2 


u 2 " v 2 


(2-65b) 


(2-65c) 


where unity scale factors are being assumed. Combining equations 
(2-65) yields 


V 1 + 2 *21 V 1 v 2 + (1 + *21> v 2 * u o 
which Is recognized to be the equation of an ellipse. 
(2-66) in polar coordinates, let 


( 2 - 66 ) 
To express 


v^ » p cos 0 


v 2 - p sin 0 


Then 


solving for p 
2 


u 


2 2 

1 + a 21 sin 20+ a 21 sin 


To find the extreme values of p, noting that p > 0, 
l 21 


2 2 

.2 u (2 a 0 . cos 20 + 2a>. sin 0 cos 0) 
ap o zi _ 

d0 ” 


(1 + a 21 sin 20 + a 21 sin^ 0)^ 


The stationary points of (2-70) are 

0 - - \ tan” 1 (^-) + , n- 0,1, 2, 3 


(2-67a) 

(2-67b) 


p^ £cos^ 0 + 2 a^ cos 0 sin 0 + (1 + aij^) 8 * n ^0j ■ (2-68) 


(2-69) 


(2-70) 


(2-71) 
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where the principle values have been taken. Evaluating the equation 
d 2 2 

of — for the values of 0 which are solutions of equation (2-71), 
it is found that the maxima of p correspond to n even and the minima 
to n odd. 

Thus if the transform of the reference data scene is scanned by a 
constant frequency aperture, the corresponding search aperture on the 
transform of the distorted data must necessarily be elliptical. Thus 
the function to be minimized is 


N , 9 

e - S | w rJ .(p) - w fei (p )) 


i»l 


(2-72) 


where 


'*[77 


~ 2 7~T~ 

a 21 sin 20 + a 21 sin 0 


■] 


1/2 


(2-73) 


for the coefficient a 


21 


(3) Rotation and Skew 

These distortions are coupled in the sense that it is not possible 
to solve for each individually. Thus the equation that must be 
minimized is 

N , 

e - E |w ri (p,0) - w bi (p',0')p (2-74) 

i”l 


[ 


1 + a 21 sin 20' + a 2 j sin 2 0 


r] 


1/2 


0 ' “ 0 + 0 

o 


where 



This is a two parameter search in which the loci are ellipses and the 
initial starting points are displaced by 0 q . As the data is experi- 
mental there always exists the possibility of local minima which will 
prevent the algorithm from reaching the true regional minima. 
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CHAPTER 3 

CORRELATION PROCESSORS 


3.0 Introduction 

In this chapter equations are developed for the optimum spatially 

invariant filter for the detection of a two-dimensional discrete data 

set S^ embedded in the background data S^. This filter, or processor, 

is linear although it may contain parameters which are determined by 

the solution of a set of nonlinear equations. These processor equations 

are interpreted as a crosscorrelation of the background data set S^ with 

a reference template obtained from the solution of a system of equations 

involving the reference data covariance matrix and the reference data 

S . 

— r 

The processor is derived for both noise free data and for data 

containing uncorrelated additive noise. All data are assumed 

distortionless since the geometrical distortion is to be removed by 

the previous data pre-processing algorithm. 

Following the derivation of the equations for this optimum filter 

several algorithms for implementing the processor are developed. These 

algorithms reflect the varying assumptions regarding the form of the 

covariance matrix of both the reference data S and of the noise 

r 

process N. The chapter concludes with a discussion of the processing 
gain given by the correlator processor. 
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3.1 Noise Model for Multi temporal Data 

A model for the noise components in the remotely sensed data is 

shown in Fig. 3-1. It is assumed that the noise, which is taken to 

be the differences between the reference data set S and its corres- 

— r 

ponding subimage in the background data S^, can be modeled as an 
uncorrelated additive source N. This source consists of two components; 
the first, N(x), is postulated to characterize such effects as sensor 
noise, non- zero background temperature and video signal quantization 
noise. The second term, D(At, X, x), is assumed to represent the 
changes in the scene due to natural processes such as the growth and 
maturing of vegetation. This latter term is of little significance 
in areas typified by the desert and mountainous regions of the south- 
western United States. However, for areas engaged in agricultural 
activities, the changes Induced by this term can be significant over 
a period of a few weeks. Thus the D component contains the information 
made available by the multitemporal data to the scientists concerned 
with the ground cover changes as a function of both the time between 
measurements and the period within the growing season that these 
measurements were taken. 

When amplitude correlation is used to measure similarity for 
the detection of two corresponding subregions of data, these temporally 
Induced changes in the scene decrease the correlation between sub images 
which are known to be geometrically corresponding. In this study the 
time between observations will be limited such that it is reasonable 
to assume that the two data sets remain partially coherent. This loss 
of correlation, however, leads to the requirement of processing a 



Physical 


Process 


r 


N - N(x) + D(At , X , x) 


Figure 3-1 A Model of Noise Processes 
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larger set of data to maintain a specified average correlator output 
signal-to-noise ratio. That these assumptions are Indeed reasonable 
rests with the experimental results presented in the sequel. 

3.2 Processor Development 

In detection theory literature the binary detection problem is 
most often developed in terms of simple hypothesis testing; the likeli- 
hood ratio is selected as the sufficient statistic and some criterion 
of optimal performance, typically the Neyman-Pearson or Bayes 
criterion, is stated. Under the assumption of additive, uncorrelated 
Gaussian noise, the resulting analytical expressions can be interpreted 
as a correlation processor. 

For the problem considered in this investigation, however, it is 
not straightforward to state the two alternative hypotheses. Rather 
the approach taken in this study is to maximize a performance index 
subject to certain constraints. The resulting analytical expressions 
again are of a form to admit a correlation implementation. 

The essential idea embodied in this performance index is the 
desire to process the background data such that the processor out- 
put corresponding to the correct superposition of j> r on _S' c is 
maximally dlscriminable from the output at all other spatial juxta- 
positions of the reference template on the background data. The 
desired processor output, therefore, consists of a large central peak, 
corresponding to the correct superposition of the template, surrounded 
by a non-zero pedestal. The filter output in the pedestal region is 
due to the interaction of the template with the background data for 
Juxtapositions of the template spatially separated from the location 
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of the central peak by distances greater than the reciprocal of the 
spatial frequency bandwidth of the scene. Such a desirable output 
Is shown In Fig. 3-2. 

The discrimination between peak and pedestal Is made by a 
comparison of amplitudes. Thus for maximum discrimination, assuming 
this central peak to be normalized to unit value, the peak amplitude 
of the filter output in the pedestal region must be minimized. Peak 
values are In general not viseful for analytical purposes since such 
a value may not having meaning. 

The peak values of the processor output can be probabilistically 
related to the variance of the process. If the output in the pedestal 
region is modeled as a Gaussian process, then the fraction of the time 
that a specified output value ia exceeded is the area in the tail of 
the Gaussian probability density function; thus, if p(x) * N(0,1), 

Prob [X > x] ■ | erfc(x) (3-1) 

where erfc(x) is the complementary error function. 

If the Gaussian assumption cannot be made but the random 
variable X has a finite second moment, then the Chebychev inequality 
may be used. Thus if 

E[X 2 ] < - , (3-2) 

then 

Prob [lx- E[X1| > c ] 1 Y ar . ^ (3-3) 


for any c > 0. 
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The thrust of this argument is that if a linear processor is 
derived such that the variance in the pedestal region is minimized, 
then it follows that the fractional time that the output exceeds the 
specified value x - x q , will also be minimized. 

To obtain an expression for this optimal processor, the quantity 
to be maximized is the ratio of the output at the point of correct 
superposition of 5» r on to the variance at all other juxtapositions. 

In the actual implementation of this filter, this operation corresponds 
to maximizing the ratio of the peak of the correlation surface to the 
variance of this surface in the pedestal region. The correlation surface 
is obtained by the crosscorrelation of the background data with the 
template generated from the reference data S^_. 

3.3 Noise Free Processor 

For noise free data the reference data is a subregion of the 

background data S^. It is assumed that the background data are a 

discrete set obtained from the product of a sample function from a real, 

wavenumber limited, two-dimensional random process with the finite 
2 

impulse grid tt defined by equation (2-2). The statistical properties 
of the reference data set are assumed to be characterized by the 
covariance matrix K. 

Using the linear spatially invariant filter H is derived such 
that when the background data are convolved with H 

z(u,v) - H * S. (3-4a) 

“0 

“ E I 8<U - i + 1, v - j + 1) h(i, j) 

i J 


(3-4b) 
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the ratio 

2, n 

z (u ,v ) 

^ ■ y« fz f u. v )] «' 5 > 

is maximized. The quantity z(u f , v ) is the filter output when the 
template is positioned on the reference data S^. This filter is 
further subject to the following constraints. The peak output 
z(u f , v ) is constrained to have the maximum value of unity. Thus 

X 1 " z * u r’ v r^ “ 1 * (3-6a) 

and this constraint can be expressed alternatively as the inner product 
relation 

“ Q!.S r ) - 1 (3-6b) 

The filter must have finite energy to be realizable, thus 

I 2 - I Z | h | 2 (3-7a) 

i j J 

- (H, H) < M < ®. (3-7b) 

3.3.1 Filter Equations for the Noise Free Processor 

The background data set ■ {s^j) * 8 a8Sume d to be an M x N 
digital picture function, and the detection template H is a P x Q 
matrix of real numbers. For this problem to have physical meaning 
the matrices must have the relationship P < M, Q < N. The filter 
output is given by the two-dimensional convolution of this picture 
function and the template 

z(u,v) ■ Z Z s(u -i+1, v-j+1) h(ij) 

i j 


(3-8) 
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where 

1 < i < P 

1 < j < Q 

and 

P + 1 < u < M + 1 

Q + 1 < v < N + 1 

These indices are chosen so Chat the two matrices are fully juxtaposed, 
thus deleting the need for consideration of end effects. This convolu- 
tion is illustrated in Fig. 3-3. 

Equation (3-8) is expressible in terms of matrix operations which 
allows the convolution to be written in terms of vector quantities; 
note there are 

m m (M - P + 1) (N - Q + 1) (3-9) 

distinct juxtapositions of the template on the background scene; each 

til 

of these superpositions is called an observation. The t observation 
is the subimage of the background data defined by the set 

ts t } - £ s t j : A < i < P + A - 1, k<j<Q + k- 1) (3-10) 

with 

1<£<M-P + 1 
l<k<N-Q+l 

The set £'s ) can be written as the equivalent pxl column vector 



Figure 3-3 Convolution of the Background Data Set 
with the Detection Template 
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u t ) - i t 




8 (t) 
”T> 


(3-11) 


where each element of si t is a q x 1 column vector 


s t (t) 


s u (t) 


L“iq (t> 


(3-12) 


The data subimage which constitutes an observation and the organization 

of the elements of this observation into a vector are shown in Fig. 3-4 

From the observations defined by equation (3-10), the pq x m 

matrix S is formed 
— a 


S 

—a 


SpU) 


Vi (1) 


s^l) 


8 ,( 2 ) 

~P 


8 (m) 
HP 


-1 


(3-13) 


tli th 

where the i column is composed of the elements of the i observa- 
tion and m is the total number of observations. 

Let 


H 


, h , 
L ~P J 


(3-14) 
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First Element in 



Figure 3-4 Organization of an Observation into an Equivalent Vector 
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be a multiple channel filter with p-input channels and one output 
channel. The filter coefficients for each channel h^ are given 
by the q x 1 column vector 

! V 1 '] 


h q d) 


The filter output z(u) can be expressed as the convolution of 
the vector quantities 


z(u) » £ js _ i+1 ( u ) u * l,...,m 


Forming the m x 1 column vector 


the convolution can be written in the matrix form as 


T 

Z - S H 
— —a — 


The ratio that is to be maximized is 


Var ^ E[(S T h) t (s t H)] 


Z Z M u z 2 (r) 

e[h t s s t h] h t k h 


(3-19) 
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T 

where K * E[S £> ] is the covariance matrix of the process and z(r) 
is the output corresponding to the correct superposition of the template 
and reference data set S^. 

To maximize this expression, z(r) will be constrained to have a 

T 

value of unity and the quadratic form H KH will be minimized. 
Additionally, from equation (3-7), the resulting H is constrained to 
have finite energy. Thus the functional relating these quantities is 
written as [25, 65] 

T i| T . 

1 “ H K H + H + X 2 HI, X 2 > 0 (3-20) 


where X^ and X 2 are the Lagrange multipliers. The extremum of this 

functional is found by first finding the stationary point of I with 

respect to H. Then noting that I is convex, it follows that this 

stationary point is the minimum of I [47]. Thus 

||-2KH + X 1 ^. + 2X2 5-0 (3-21) 


Solving for H 


H - - \ X X (K + X 2 I)" 1 S r 


(3-22) 


The values of the Lagrange multipliers X^ and X 2 are obtained by 
solving the simultaneous set of nonlinear algebraic equations 





H - 


(3-23a) 


M 


(3-23b) 
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Substituting equation (3-22) into (3-23), yields 

" 1 X 1 -r ( - + X 2 - )_1 £r " 1 (3- 24a) 

T 

| xj s£ [(K + X 2 I) -1 ] (K + \ 2 I)" 1 S r - M, X 2 > 0. (3-24b) 

The solution of the equations (3-22) and (3-24) then yields the 
solution for the desired processor. 

3.3.2 Some Observations Regarding the Noise Free Processor 

Two salient points should be noted from the form of equation (3-21). 

First, the covariance matrix is the same for both the reference data S 

— r 

and the background data since it is assumed that the process is 
stationary. Since there are no energy differences between the reference 
and background such as required for incoherent detection, this filter 
must be a coherent processor. 

Secondly, there is a striking similarity between the filter given 
by equation (3-21) and the filter for detecting stationary targets in 
a clutter background in the radar problem [61], In particular, for 
\ 2 - 0, the radar filter is the "inverse filter" which was first 
discussed by Urkowitz [59], For \ 2 4 0, the corresponding filter is 
the "clutter rejection filter" [18, 6l] which includes the considera- 
tion of additive receiver noise. The effect of these filters in the 
radar problem is to increase the bandwidth of the signal being processed; 
it has been shown that for the constraint of a fixed signal bandwidth, 
the signal spectrum for the optimum detection of targets in a clutter 
background is flat over this specified bandwidth [18]. 
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The one-dimensional analog of the filter given by equation (3-21) 
for X 2 “ 0 is also well known in the field of seismic data processing. 

In this area the filter is known as the "spiking" filter. 

3.4 Noise Free Algorithm Development 

In this section use is made of both equation (3-22) and various 
assumptions regarding the form of the covariance matrix K to develop 
algorithms for implementing the processor on a digital computer. Care 
must be exercised in discussing the stationarity of the process 
characterized by the covariance matrix. In one- dimensional data, a 
result of stationarity is the requirement that the covariance depends 
only on the separation of the points t « - x^ However, in two 

dimensions this definition is too restrictive in that it is conceivable 
that the covariance can depend on both the spatial separation t ■ x^ " 
and the direction of the line passing through these points. 

Algorithms are developed for both cases; however, for non- isotropic 
covariance matrices a general covariance matrix must be considered. 

This in turn requires the solution of a large set of linear equations 
for determining the template, thus limiting the value of this procedure 
for reasons of computational tractability. 

Major emphasis is given to algorithms for X 2 “ 0. Recall from 
equation (3-7) that this constraint was related to the filter energy. 

For finite dimensional filters this energy is finite so long as the 
elements of H are finite; this must be true if the covariance matrix 
is non-singular. In addition for the noise free processor, the 
variance in the output is due to the Interaction of the reference 
template with the background signal. Thus the correlator output 
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signal -to-noise ratio, as defined by equation (3-19), is invariant 
with respect to the filter energy. 

The algorithm for the general covariance matrix is considered 
first and, following this, study will be given to more restrictive 
assumptions regarding the form of the covariance matrix. In this 
latter case, the matrix assumes the form of a block Toplitz matrix 
which has much practical interest. 

3.4.1 General Covariance Matrices 

If the two-dimensional process from which is selected has 

the general PQ x PQ covariance matrix K, then the solution of the 
processor template H is given by equation (3-22) 

H - - ■§ *! K" 1 S r (3-25) 


and from equation (3- 6b) 


1 T -1 
-|Xi S*K 1 S r - 1 


(3-26a) 


and solving for X. 


-2 


1 T -1 
S K S 


(3-26b) 


Thus the template for the crosscorrelator is given by 


K _1 S 


H 


T -1 
S K S 

— r — — r 


(3 -27a) 


®1 S' 1 — r 


(3-27b) 
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where 

g x - <4 S’ 1 s r ) _1 (3-28) 

is Che normalization coefficient. 

To estimate the covariance matrix K, the data in the region 
surrounding S is partitioned into the subregions as shown in Fig. 3-5. 
Each subregion is called an observation and a correspondence is set up 
between the elements of an observation and the vector or as 

•ij 4* <^(1) (3-29) 

as illustrated in Fig. 3-5(b). Thus the i^ observation is mapped into 
the vector a (i). 




(3-30) 


The covariance matrix is then estimated by the unbiased estimator 


" M 

E a (i) a (i) 
i-1 



M 

E a (i) a (i) 
, n 1 


M 

I <*i(i) <*(i) 
i-1 


M 

E a (i) a (i) 


(3-31) 
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(a) Partition of Data Region into Observations 
for Covariance Estimation 



(b) Packing Data Observations into a Vector 


Figure 3-5 Data Format for Computing Covariance 
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where M is the total number of observations and n is the number of 
elements in a. If, for example, the reference data set is a subimage 
of 32 x 32 elements, then the resulting covariance matrix is 
1024 x 1024 elements in size. 

Premultiplying equation (3-27b) by K yields 

KH-g^ (3-32) 

Reference data array sizes generally range between 32 x 32 to 
128 x 128 elements in practical algorithms using experimental data. 

Thus to solve equation (3-32) for data arrays in this size range requires 
a computer with the capability of handling very large array sizes. 


3.4.2 Covariance Matrices of the Toplitz Form 

The stationary isotropic covariance matrices represent an interest- 
ing class of data for filter development; the symmetry properties of 
these covariance matrices can be utilized to reduce the complexity of 
the processor. Covariance matrices for a process of this class have 
the form 




(3-33) 


L K. # • • h. 

— m -o J 

This form of matrix is called a block Toplitz matrix and its properties 
have been investigated by Grenander and Szego [3l], An algorithm for 
inverting a matrix of this form was first given by Levinson [45] and 
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such inversion algorithms have more recently been studied by 

Robinson [55, 62] and by Bednar and Farmer [8]. The principal 

advantage of this algorithm is that the number of operations required 

2 3 

for inversion is proportional to m rather than m as in conventional 
me thods . 

For each observation of the data arranged as shown in Fig. 3.5(b), 
the covariance matrix has the partitioned form 
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where a unit variance process is assumed. It is immediately evident 
that the matrix (3-34) has the same form as matrix (3-33) with each 
partitioned block of (3-34) corresponding to the data obtained from 
one line of the observation shown in Fig. 3-5(b). 

With experimental data it is unreasonable to assume that the 
blocks along a diagonal will be equal. However, if the process is 
assumed to be of this class, then the blocks are averaged to provide 
a matrix of the required form. The template is again obtained from a 
solution of equation (3-32) , 

The difficulty in applying this algorithm is that a large set of 
data is required for estimating the covariance matrix. However, it 
has been observed that the correlation function for a large number of 
observations of remotely sensed data exhibits an approximate exponential 
form [52]. Thus, it is not clear that this algorithm would yield 
superior results to the next algorithm to be developed. 


3.4.3 Isotropic Exponential Covariance Matrices 

A very simple processor results if it is assumed that the statis- 
tical properties of the reference data set are modeled by an 
isotropic exponential covariance matrix. This matrix has the block 
Toplitz form with the diagonal terms given by 

k u - || PlJ li (3-35) 



with 


i - J 
i * 1 


(3-36) 
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where p is equal to the average adjacent cell correlation coefficient 
in the reference scene. The off-diagonal terms are given by 


k - p k 

-ij p -ii 


(3-37) 


with p defined as for equation (3-35) . 

This covariance matrix can be written as 


K - 


in 

“in 

P 2 iu 

• • • 

“°iii 

Pin 

in 

“ in 

• • • 

P- I iu 


A u 


-ii 


(3-38) 


The inverse of this matrix has been shown to be [6] 


K 


1-P 


hi 


'“in 


, -l 

’“in 


2 -1 
(1 + P >k 1 


ii 


'“ill 


... 0 


-“in 

(i + P 2 )k'J 


* -ii 


(3-39) 


where 


-1 1 


-ii _ 2 

1-P 


1 -p 0 

. . 2 

•0 1 + P -P 


0 -p 


1 + P 


.. 0 


... 1 


(3-40) 
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To obtain a solution for the desired template, substitute 
equation (3-39) Into equation (3-28) and neglecting the multiplicative 
scalars 


k" 1 

1*11 


-*U 


a (r) 


H - 


-pk'J (1 + q 2 )k^ -pkjJ 

0 -pk [I (1 + p 2 )k‘J 


( r ) 




(3-41) 


where the reference data set has been written as an equivalent 
column vector as defined by equation (3-11). Expanding equation 
(3-41) , the l** 1 filter element h^ Is 

h - [-*;! ii+1 + < l + P 2 )iu ii - p iii -Vi] <3 - 42) 

where the identifier that the s^ are from the reference data has been 
omitted. Now substituting equation (3-40) into equation (3-42), again 
omitting multiplicative constants, yields the elements of the template. 
Thus 

h ik - - P [-P s 1+ljk+1 + a + P 2 > » 1+l>k - P » 1+1 , k .!] 

+ (1 + p 2 ) [-P. l>k+1 + U + P 2 ). i>k - P 
- P [-P » 1 . ljk+1 + U + P ) - P s i-l, k -l] <3 ' 43) 
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Expanding, 

h ik " p 8 i+l,k+l " p(1 + p * S i+l,k + p 8 i+l ,k-l 

- »<* + p2 > -i.wi + < 1+ p2 > 2 *i.k - P<1 + p2 > •i.k-i 

+ p2 - p(1 + p2) Vl,k + p2 “l-l.k-1 

(3-44) 

Equation (3-44) is recognized as one term of the discrete convolution 
of the grid operator G, where 

T 2 , 2 V 2 

p -p(l + p ) P 

G- -p(l + p 2 ) (1 + P 2 ) 2 -P(l + P 2 ) (3-45) 

2 2 2 

L o -Pd + p z ) p J 


with the reference data set S . This convolution is written as 

— r 



" 

- — 

— 

— . 

— ' 


— 

8 i-l,k-l 

H 

1 

•H 

<0 

8 i-l,k+l 

— 

H - G * 

— 

8 i,k-l 

8 i,k 

8 i,k+l 

— (3-46) 



8 i+l,k-l 

8 i+l,k 

S i+l,k+l 

— 


where 

2 < i < P - 1 
2 < k < Q - 1 

and only the elements of S used for evaluating h . have been explicitly 

— r ix 

Indicated. The indices have been chosen to eliminate the end effects of 
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the convolution; therefore, the reference data set should be two rows 
and two columns larger than the desired template size. 

The result of the development of G is that rather than having to 
evaluate H using equation (3-41), U can be evaluated directly using 
equation (3-46). The resulting simplification is very substantial. 

The grid operator G has several interesting properties. For p - 0, 
equation (3-46) becomes 



(3-47) 


Thus if the reference scene consists of a set of uncorrelated points G 

Ss 

becomes the identity grid operator for the convolution of equation (3-46) 
and h tk « 8 lk . Thus the template is the reference data set itself. 

For p ■ 1, the grid operator becomes 



(3-48) 


This is the discrete approximation to the mixed fourth partial 
derivative, obtained by convolving the discrete approximations to 
the second partial derivative along each coordinate axis. Thus 


G - G x * G 2 - 


0 

0 

0 


0 

1 

0 

1 

-2 

1 

* 

0 

-2 

0 

_ 0 

0 

0 


0 

1 

0 


(3-49) 


The coefficient p required for determining the grid operator G 
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is experimentally obtained from the data. A useful procedure is to 
re- format the reference data S^, which is represented as the P x Q 
matrix, into the [L x 2] matrix a {Sp(i,j)] with 

L - [ f ] Q (3-50) 

where [•] is read "the greatest integer less than or equal to." 

The correlation coefficient p is then given by 


P 


L 

I » (1,1) 8 (i,2) 

-i P P 




L 

L 

i-1 





(3-51) 


3.4.4 Non-zero Constraints on Filter Energy 

The solution for the crosscorrelation template was given by 
equation (3-22) as 

H - - \ Xj (K + X 2 I)" 1 S r (3-52) 

Recall that the Lagrange multipliers X^ and X 2 were introduced into the 
problem formulation so as to constrain the filter output amplitude and 
energy respectively. To solve for X^ and X 2 it was found necessary to 
solve the simultaneous system of equations 

- \ (K + X 2 I) -1 S r - 1 (3-53) 


[(K + X 2 I)' 1 ] 1 (K + X 2 I)" 1 S r - M 


(3-54) 
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From equation (3-53) solve for X^ 


X, - -= — = (3-55) 

S (K + X 2 I) S* 


and substituting equation (3-55) into equation (3-54) yields 


T 

s l [(K + X 2 I)" 1 ] (K + X 2 I)" 1 S r 

S T (K + X. I)' 1 S 
-r — 2 — — r 


M 


(3-56) 


and further 

T 

S* [CK + X 2 I)" 1 ] (K + X 2 I)' 1 S r - M (K + X 2 I)‘ l - 0 

(3-57) 

This is a nonlinear equation in the one known X 2 - To solve this 
equation requires an algorithm for evaluating (3-57) which includes 
a nonlinear root finding program. If the value of X 2 obtained from 
the solution of equation (3-57) is negative, the constraint is not 
active for the specified value of M and X 2 is set equal to zero [65]. 

3.4.5 Signal Prewhitening 

The function of the Inverse covariance matrix is to prewhiten both 

the input scene and the stored reference picture function. A processor 

for the detection of in the background ^ is shown in Fig. 3-6(a). 

In this implementation the observations £ from are crosscorrelated 

with W S where 
— r 

w s - k ’ 1 s 

— r — — r 


(3-58) 









87 


An equivalent expression is obtained if K * is factored into the upper 
triangular matrix and its transpose 

K" 1 - W W T - W T W (3-59) 

which leads to the processor illustrated in Fig. 3-6(b) . 

It should be noted that the elements of W S of Fig. 3-6 (b) are 
uncorrelated and are equivalent to samples from a unit variance white 

V 

noise process. This is shown by evaluating the covariance of £> r ; 
thus 

E [ W S r W Sy ]- ^ (W sf] - w E[s r ] W T 

- W K W T - W( W T W) “ 1 W T - I (3-60) 

The prewhitening of the observations of is illustrated in 
Fig. 3-7. The scatter diagram exhibiting the correlation between 
adjacent elements of is shown in part (a) of this figure. It is 
Inferred from equation (3-60) that an equivalent representation of W S f 
is a vector consisting of samples from a unit variance process; the 
scatter diagram of adjacent elements of this processed vector is 
illustrated in Fig. 3-7(b). To be noted is that the variance along 
each coordinate axis is equal, resulting in a circular scatter diagram. 

When the data become highly correlated the points in the scatter 
diagram tend to become concentrated along the unity correlation axis; 
thus the restructuring of the data into the circular scatter diagram 
requires a gain coefficient greater than unity along the axis per- 
pendicular to the unity correlation axis. For values of data 
correlation approaching unity, it is reasonable to assume that the 
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algorithm will introduce serious numerical error into the processed 
vector and this resulting vector will be of little value. 

Assuming that the covariance matrix of a data set is known, it 
is suggested that a measure suitable for determining if a vector from 
this data is amenable to prewhitening is the determinant of K, or 
equivalently det(K ^). For a unit variance process the value of 
det(K *) must always satisfy the relation. 

det(K _1 ) > 1 (3-61) 

To show this, first note that 

det (KK' 1 ) - det (I) -det (K) det (K _l ) - 1 (3-62) 

thus, 

<s _1 > - < 3 -«> 

A covariance matrix has the properties of being symmetric and positive 
definite. Thus the matrix K is diagonalizable by the similarity 
transformation [26]. 

diag (K) - P _1 K P (3-64) 

where P is the matrix whose columns are the eigenvectors of K and 
the resulting diagonal elements of diag(K) are the eigenvalues X^ of K. 

This matrix diagonalization is the well known principal component 
transformation [36], and in general the eigenvalues can be ordered 

- * * * - X n > 0 (3-65) 

The energy in the process is represented by the sum of the eigenvalues 
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n 

E - £ X (3-66) 

i-1 

The determinant of the diagonal matrix is 

n 

det [diag (K)] - n X (3-67) 

i-1 

and by the result derived in Appendix 0 

det [diag (K)] < 1 (3-68) 

which implies that 

det [K] < 1 (3-69) 

since the value of the determinant of the matrix K is invariant under 
the similarity transformation P * K P. 

Thus, it is concluded that 

det (K" 1 ) > 1 (3-70) 

The equality holds only for the identity matrix which represents an 
uncorrelated process. As p -* 1 however, the value of det(K *) becomes 
large because at p - 1, K becomes singular and det(K) - 0. 

From these considerations it is concluded that there must exist 
some range of det(K *) such that 

1 < det(K _1 ) < M x < « (3-71) 

where the algorithm will work effectively. For those data sets where 
det(K ■*■) > Mj, the data become too highly correlated to be of value and 
consequently should be rejected as a data set to be processed for 
obtaining registration. For the isotropic data sets this becomes a 
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matter only of rejecting those data for which the average adjacent cell 

correlation p is greater than some p 

6 K max 

3.5 Multitemporal Data Processor 

In processing multitemporal data there must be available a 

minimum of two data sets; as these data sets are processed pairwise, 

the discussion will be phrased in terms of two sets. Subsets of data 

are selected from one set and are used as the processor reference 

signals. The background data are chosen from the second data set. 

It is assumed that the reference signal S is noise free. It is 

— r 

further postulated that the approximate location of S^, the corres- 
ponding data set in the background data set, is known, and that the 
differences in amplitude between S^ and s/ are assumed to be due to 
additive noise These assumptions allow the development of the signal 
processor in an analogous manner to that done for the noise free case; 
the major difference is the addition of a covariance term reflecting 
the addition of the noise. 

3.5.1 Filter Equations for the Multi temporal Processor 

In a procedure analogous to that employed in section 3.3.1, the 
constrained problem is formed with 

S - S + N (3-72) 

— s — r — 

The expected value of the output is constrained as 

I r - E[z(r)] - E[H T (S r + N ) 1 “ H T S r = 1 (3 _ 73) 

with E[N] = 0. The filter is to have finite energy so that 
! 2 = H T H < M < ® 


(3-74) 
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The objective is to maximize the ratio given by equation (3-19) by 
minimizing the variance of the filter output 


Var f Z] - E 

T 

[h (S + N)] H (S + N)| 

(3-75a) 

T 

■ H 

E [ (S T + N T ) (S + N)] H 

(3-75b) 

T 

- H 

(K + K ) H 
—8 — n — 

(3-75c) 

E[S T S] - K 
m " 1 8 


(3-75d) 

E[N X N] - K 

— — — n 


(3-75e) 


E[S T N] « E[N T S] =0 (3-75f) 

Thus all of the equations of subsection 3.3.1 are now modified by 

\ 2 > 0 (3r76) 

(3- 7 7a) 

(3- 7 7b) 
(3-77c) 

As noted previously in the development of equation (3-21) , the 
functional I is convex. Thus the solution for H, equation (3-78), gives 
the stationary point of 1 and this solution is the minimum. 


replacing K by K + K . Thus 

— — s — n 

t - h t (k + K )H + x. h t s + X. h t h, 

— —8 — 1 — — r 2 

which has the solution for the stationary point 

2 * - i h <*. + £n + X 2 !>'* — r 

with 

T 

H S - 1 

— — r 

H T H < M 
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3.5.2 Multi temporal Algorithm Development 

For the multitemporal processor both the reference and noise 

covariance matrices must be determined. Assuming there is sufficient 

data for computing these matrices, the reference covariance matrix, K , 

1 8 

is computed from the reference data as suggested in subsection 3.4.1. 

The additive noise component is taken to be the differences between the 
reference data S f and the corresponding data set in the background 
data. The equation of this noise component is written 

N - S' - (1 - k) S r (3-78) 

where a differential gain k has been considered. To determine the 
coefficient k, the stationary point of 

e 2 - N T N - [s; - (1 - k) s/ [S' - (1 - k) S r ] (3-79) 

with respect to k is found. Thus 
2 

- 2 [S^ - (1 - k) S/ S r - 0 (3-80) 

which has the solution 

-T 

S S 

k - 1 - —r 11 (3-81) 

S S 
— r — r 

2 

This value of k yields the minimum of e since this functional is 
convex. 

The noise covariance matrix K can now be computed. With Ki»K + K 
the algorithms of section 3.4 are applicable to the multi temporal 
processor. 
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The difficulties encountered in the solution of equation (3-77) 

for the multi temporal processor include those requirements observed in 

subsection 3.4.1, the necessity for handling very large arrays. 

However, a more fundamental problem is the unavailability of a 

sufficiently large set of data, having similar statistical properties 

as the reference data, such that a useful estimate of K can be made. 

— s 

An additional difficulty is encountered in estimating the noise 
properties. With the noise defined as the differences between and 
s/, the determination of requires that these data be essentially in 
registration such that the differences are due to scene changes and 
are not the result of misregistration error. 



(3-86) 
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For the noise free case. K - 0 

— n — 

I “ 4 C *r < 3 ' 87a > 

- Tr (S K -1 S T ) - N (3-87b) 

— r —s ~t 

where N is the number of observations averaged. For K 4 0 

— n — 

I - 4 <£„ + y 1 i r (3-88) 

In the derivation of an equivalent result by Arcese et al. [6], 
there is an error in the derivation. A development such as suggested 
by these authors is given in Appendix C and the result is shown to be 
consistent with equation (3-87). 
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CHAPTER 4 

EXPERIMENTAL RESULTS 

4.0 Introduction 

In this chapter certain of the algorithms developed In the 

previous chapters are applied to both Ideal data sets and to Imagery 

obtained with a line-scanning spectrophotometer mounted in an aircraft 

The algorithms used with a specific data set, and taken collectively, 

are termed a processor. The results obtained using the ideal data 

are presented first, followed by the application of the processor to 

a scanner data set in which the reference data S is selected as a 

r 

subimage of the background data. This, of course, yields a distortion 
less and noise free reference data set. 

Multitemporal data sets are then considered. Search techniques 
utilizing the modulus of the Fourier transform of these data are used 
to estimate the regional geometric distortion, and results of the 
processor operating with these corrected data are obtained. 

4.1 Ideal Data 

The purpose of generating ideal data sets is to obtain picture 
functions with known statistical properties. These data are then used 
to both check for the correct operation of the computer codes as well 
as to evaluate the performance of the processing algorithms with 
variations in the data parameters. These ideal data were obtained by 
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convolving a two-dimensional array of computer generated random 
numbers, having a Gaussian distribution, with a two-dimensional low- 
pass filter having a specified impulse response. 

Two sets of ideal data were generated, differing in the form 
of the filter used to smooth the set of random numbers. The impulse 
response of the first of these filters was chosen so as to obtain 
a data set with a covariance function modeled by equation (3-38); 
this data set exhibits a strong correlation of the data along the 
coordinate axes. An alternative data set having no preferential 
direction of data correlation was obtained with the second filter. 

4.1.1 Separable Exponential Filter 

The impulse response of the filter used in this subsection, 
termed the "separable exponential filter," is given by 

h(x) > k 1 exp [^- (x L + x 2 )] (4-1) 

o 

|xj < R, i - 1,2 

The normalization coefficient , derived in Appendix B, is 

k l ■ { 4 R o t 1 ■ exp (r-)]} (4 - 2) 

o 

where R q is defined as the "characteristic length" of the filter. The 
variable R is the size of the truncated impulse response with both R 
and R q being measured in units of picture elements. A plot of the 
impulse of this filter is shown in Fig. 4-l(a). 

This filter was chosen so as to yield a data set with statistical 
properties closely approximating the isotropic exponential covariance 
matrix defined by equation (3-38). Thus the processor for this data 



Figure 1+-1 Correlator Outputs for Various Amounts of Data Correlation 
(Separable Exponential Filter) 





Figure 4-1 Cont . 
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set is the crosscorrelation of the background data with the template 
given by equation (3-46). This template is obtained by simply by 
convolving the grid operator G 


G - 


2 

P 

-P(l + P 2 ) 


2 

P 


2 2 
-p(l + p ) p 

2 2 2 
(1 + pV -p(l + p*) 

-p(l + p p 


(4-3) 


with the subimage of the background data selected as the reference 
data set S^,. The adjacent cell correlation coefficient p of the 
reference data J> r is used to determine the grid operator, and 
experimentally observed values of p for various values of R q and 
template sizes are given in Table 4-l(a). 

The processor outputs are shown in Figs. 4-l(b)-(f) for data 
sets having specified values of R q and where the template is a subimage 
of the background data. The correlator template size for these 
results was 32 x 32 picture elements. There are several properties 
of these outputs which are of interest. The most obvious result is 
that the output function becomes smoother and the pedestal variance, 
relative to the central peak, Increases as R q increases; an increase 
in R q decreases the spatial frequency bandwidth of the two-dimensional 
signal. Thus as the spatial -frequency bandwidth is reduced, there is 
both a loss in signal- to-noise ratio, as defined by equation (3-19), 
and a loss of resolution. In this application resolution is defined 
as the width of the central correlation peak at one-half its height. 

It is also observed that the peak values of the output in the 
pedestal region appear to be distributed such that there are no large 
distinct peaks in this region. Thus, a unique central peak can be 
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Table 4-1. Experimental Values of the Average Adjacent 
Cell Correlation Coefficient (p). 

(a) Separable Exponential Filtered Data 


Template Filter Characteristic Length (R q ) 


Size 

0.0 

0.5 

1.0 

2.0 

4.0 

32 x 32 

.049 

.309 

.652 

.868 

.953 

48 x 48 

-.009 

.296 

.647 

.874 

.954 

64 x 64 

-.013 

.297 

.650 

.887 

.962 

(b) Isotropic Exponential Filtered Data 



Template 

Filter Characteristic Length 

<V 

Size 

0.0 

0.5 

1.0 

2.0 

4.0 

32 x 32 

.049 

.254 

.658 

.902 

.972 

48 x 48 

-.009 

.260 

.676 

.904 

.964 

64 x 64 

-.013 

.268 

.699 

.902 

.967 
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unambiguously located. 

Ideally the processor for the noise free signal given by equation 
(3-27) is the inverse of the filter used to obtain the data. Defining 
a parameter Of with 

0 < a < 1 (4-4) 

as the prewhitening coefficient and choosing 

p y = or p (4-5) 

as the parameter for the grid operator, equation (4-3), yields the set 
of outputs shown in Fig. 4-2. The background data set used for this 
sequence of correlator outputs was obtained with a filter having the 
parameter R q ■ 2, and the correlator template was 32 x 32 picture 
elements in size. 

For the case of or - 0.0, the processor reduces to a cross- 
correlation between the background data and a template consisting of 
a subscene of S_^. For a - 1.0, the processor ideally implements the 
inverse filter. The effect of the prewhitening coefficient for or < 0.5 
is observed to be small. However for or > 0.7, the processor template 
begins to significantly alter the correlator output. 

Experimentally obtained curves giving the output signal- to-noise 
ratio as a function of the prewhitening coefficient are given in 
Fig. 4-3. This family of curves for various values of filter character- 
istic length was obtained by computing the processor output for each of 
the set of data parameters (R q , or) where the processor template was 
32 x 32 picture elements in size. The signal-to-noise ratio of each 
of these correlator outputs was then determined. The signal component 



Figure 4-2 Correlator Outputs for Various Prewhitening Coefficients 
(Separable Exponential Filter) 






Figure 4-2 Cont. 





Output Signal-to-Noise Ratio 
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Alpha (a) 


Figure 4-3 Correlator Output Signal-to-Noise Ratio as a 
Function of the Prewhitening Coefficient 
(Separable Exponential Filter) 
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Is defined as the square of the difference of the amplitude of the 
central peak and the average value of the output in the pedestal 
region. Thus 

S„ - r*(£) - <*(u)> ~| 2 (4-6) 

L 

The noise component is defined as the variance of the correlator in 
the pedestal 


N 

c 


1 

N -1 
P 


N 

E 



<*<u)> ] 

ujf£ J U)lr 


(4-7) 


where is the number of points used in the summation. It should be 
noted that the condition u ^ r used in equations (4-6) and (4-7) 
implies that the central peak was deleted and not just the point u - r. 
The average value was also taken of the set of points having the 
central peak deleted. 

The significant results to be obtained from these plots are that 

the signal- to-noise ratio increased approximately 6 dB for reasonably 

correlated data and that the slope of each of the curves near a <= 1.0 

is relatively small. This latter observation means that the processor 

is relatively insensitive to the choice of p for values of p near the 

correct value P^. The performance of the processor for each value of 

R is given in Table 4-2. 
o 

The processor output signal- to-noise ratio, given by equation 
(3-88b), is theoretically independent of the adjacent cell correlation 
coefficient. However, for the class of data being considered in this 
section, as p -* 1 the amplitude of the data must approach a constant 
value. Curves giving experimental results of the output signal-to-noise 
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Table 4-2. Correlator Performance for the Separable 
Exponential Filtered Data. 


R 

o 

P 

0 

1 

O 25 

G a-1 

Processor 
Gain (dB) 

0 

.049 

1006 

1001 

0 

0.5 

.309 

773 

1001 

1.1 

1.0 

.652 

224 

878 

5.9 


2.0 


868 


52.6 


441 


9.2 
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ratio as a function of p are given in Fig. 4-4. It is observed that 
the correlator output signal- to-noise ratio decreases by approximately 
3 dB at p ■ 0.9. For values of p > 0.9 the signal-to-noise drops 
precipitously. It is concluded that the processor yields acceptable 
results for the range 0 < p < 0.9, and the value of the results for 
p >0.9 is questionable. 

4.1.2 Isotropic Exponential Filter 

The experimental procedure of subsection 4.1.1 is repeated in 
this subsection where the data are obtained by convolving the two- 
dimensional array of random numbers with a circularly symmetrical 
exponential filter. The impulse response of this filter is of the form 


h(x) - k 2 exp [- ~ (x 2 + x 2 V /2 ] 

o 

(4-8) 

|xj < R, i - 1,2, 



where the normalization coefficient derived in Appendix B, is given 
by 

k 2 ■ { R o t 1 - (1 + F> «■> <- I"*] } 1 <4 - 9 > 

o o 

The parameter R q is defined as the characteristic length of the filter 

as was done earlier in equation (4-2). The variable R is the size of 

the filter, and both R and R are measured in terms of picture elements. 

o 

A plot of the impulse response of this filter is shown in Fig. 4-5(a). 

The data set used in subsection 4.1.1 had preferential correlation 
directions introduced by the filter defined by equation (4-1). The 
filter used in this section was selected so as not to introduce such 


Correlator Output Slgnal-to- Noise Ratio 


9 



Average Adjacent Cell Correlation (p) 


Figure 4-4 Correlator Output Signal-to-Noise Ratio aa a 
Function of the Average Adjacent Cell 
Correlation 

(Separable Exponential Filter) 




(a) Filter Impulse 
Res ponse 


(c) R 

o 


Figure 



-5 Correlator Outputs for Various Amounts of Data Correlation 
(isotropic Exponential Filter) 






Figure U-5 Cont . 
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preferential directions into the correlation function; this filter is 
termed the "isotropic exponential filter." Observed values of p for 
various values of R q and template sizes are given in Table 4-l(b). 

The purpose of this experimental procedure is to compare the 
results of subsection 4.1.1 with the processor performance when the 
data set has statistical properties other than those given by 
equation (3-38). The correlator outputs for corresponding values 
of R q as used in Fig. 4-1, but using this isotropic set of data as 
the input, are shown in Figs. 4-5(b)-(e). The template size was 
again 32 x 32 picture elements. Comparing the corresponding outputs 
in Figs. 4-1 and 4-5 shows that the outputs have very similar character- 
istics. 

Defining p' and a as in the previous section, the sequence of 
processor outputs as the prewhitening coefficient a is varied, is 
given in Fig. 4-6. The data used in this sequence were obtained with 
the filter having R q - 2.0 and the correlator template size was 32 x 32 
elements . 

The numerical results of the two sequences of experiments shown 
in Figs. 4-5 and 4-6 are summarized by the curves of Fig. 4-7. The 
identical procedure for determining the signal-to-noise ratios as used 
to obtain Fig. 4-3, was used to compute the curves of Fig. 4-7. The 
performance of the processor with the isotropic data set is summarized 
in Table 4-3. 

Care must be used in comparing the results shown in Figs. 4-3 

and 4-7 because the value of the average adjacent cell correlation 

coefficient differs between the data sets for the same value of R . 

o 



Figure U-6 Correlator Outputs for Various Prewhitening Coefficients 
(isotropic Exponential Filter) 






Figure b-6 Cont. 
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Alpha (a) 


Figure 4-7 Correlator Output Sigual-to-Noise Ratio as a 
Function of the Prewhitening Coefficient 
(isotropic Exponential Filter) 
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Table 4-3. Correlator Performance for the Isotropic 
Exponential Filtered Data 


R 

o 

P 

G o-0 

G a(M4X) 

Processor 
Gain (dB) 

0 

0.049 

1006 

1005 

0 

0.5 

.254 

735 

967 

1.2 

1.0 

.654 

153 

813 

7.2 

2.0 

.902 

28.5 

558 

12.9 


4.0 


972 


9.58 


423 


16.6 
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However, It is readily observed that the curve for R q ■ 0.5 of 

Fig. 4-7 does not reach the curve for R ■ 0. In addition the curves 

o 

of Fig. 4-7 have their peak values at approximately or - 0.95. 

Fig. 4-9 gives the correlator output signal-to-noise ratio as a 
function of the average adjacent cell correlation for the family of 
curves generated by templates having N x N elements. Comparing Fig. 4-9 
to the corresponding earlier result shown in Fig. 4-4 yields insignifi- 
cant differences. Thus it is concluded that the processor is not 
critically dependent on the exact form of the covariance function of 
the data in order to give acceptable performance. 

Certain additional results were computed using the data set 
obtained with the isotropic exponential filter. The correlator output 
for various sizes of templates is shown in Fig. 4-8. The data set 
used was obtained using the filter with R q - 2.0 and the template sizes 
were N x N picture elements where the values of N were 32, 48 and 64. 

It is observed that the variance in the pedestal decreases as N 
increases; however, the width of the central peak remains essentially 
unchanged . 

In Fig. 4-10 the processor output signal-to-noise ratio is plotted 
as a function of the template size; in these plots the independent 
variable was chosen as N where the number of picture elements in the 
template is N x N. Fig. 4-10(a) gives the results where the template 
is a sub image of the background whereas Fig. 4-10(b) gives the results 
where the template is obtained by the processor for the prewhitening 
coefficient a ■ 1.0. It is of interest to note that the results of 
Fig. 4- 10(b) are a series of linear loci and are approximately parallel. 



Figure 4-8 Correlator Output for Various Template Sizes 
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Average Adjacent Cell Correlation (p) 


Figure 4-9 Correlator Output Signal-to-Noiae Ratio as a 
Function of the Average Adjacent Cell 
Correlation 

(isotropic Exponential Filter) 
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Figure U-10 Correlator Output Signal-to-Noiae Ratio as a 
Function of Template Size 
(isotropic Exponential Filter) 
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From the results of the two previous sections. It Is observed 
that the differences which exist between corresponding processor 
outputs for the two data sets are not significant. Thus it is suggested 
that this processor which is developed for the class of data sets having 
isotropic exponential covariance may be applicable to data having 
different statistical properties. An additional result of interest is 
the smoothness of the curves of signal-to-noise ratio as a function of 
the prewhltenlng coefficient for values of a near a - 1. This suggests 
that the choice of a value for p is not a critical factor for the 
processor to yield useful results. 

4.1.3 Signal-to-Noise Ratios for Geometrically Distorted Data 

If geometrical distortion exists between two corresponding sets 
of data, there results a loss in the correlator output signal-to-noise 
ratio compared to the value obtainable using undistorted data. An 
experimental verification of this statement was obtained and the results 
are shown in Figs. 4-ll(a) through (c) . These curves were computed 
using the ideal data discussed in section 4.1. 

The distorted data were obtained by rotating the set of computer 
generated random numbers with a nearest-neighbor interpolation algorithm. 
Following this rotation these data were convolved with the isotropic 
exponential filter, equation (4-8), having the appropriate choice of 
characteristic length R q . These data were then crosscorrelated with a 
similarly filtered set of undistorted data, and the resulting output 
signal-to-noise ratio was determined. 
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Rotation (Degrees) 

(a) R 0 - 0.0 


Figure 4-11 Correlator Output Signal-to-Noise Ratio as a 
Function of Rotational Misregistration 
(isotropic Exponential Filter) 
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The significant results shown by these curves Is that as the 
data becomes more highly correlated, the correlator output signal-to- 
noise ratio, although less in magnitude, is less sensitive to the 
geometrical misregistration. This of course is an intuitively satisfy- 
ing conclusion since the correlation distances become greater as the 
characteristic length of the filter increases. 

4.2 Noise Free Experimental Data 

The processor used in section 4.1 is now applied to line scanner 
data. In this section the reference data is again selected from 
the background data, giving a reference signal which is both noise 
free and distortionless. The Images processed in this section are 
typical of scenes observed in line scanner imagery. 

In the results of this section, shown in Figs. 4-12 through 4-14, 
the title "image template" refers to the use of itself as the corre- 
lator template, whereas the term "processor template" is taken to mean 
that £> r has been convolved with the grid operator G, given by equation 
(4-3), to obtain the template for the correlator. In addition the 
number of picture elements in the template for a given correlator 
output is N . 

A scene of an agricultural area is shown in Fig. 4-12(a); the 
outlined area represents the approximate area from which the reference 
data set was selected. Fig. 4-12(b) shows the output of the cross- 
correlation of with the background j>, . It is noteworthy that 
Fig. 4-12(b) typifies the correlator output of those scenes consisting 
primarily of agricultrual fields. The correlation function is being 
viewed at an angle of tt/ 4 radians, and it is readily apparent that the 



(a) Agricultural 
Area 


(c ) Processor Template 
N 2 = 4096 


(b) Image Template 

N 2 = 4096 


Figure 4-12 Correlator Output for an Agricultural Area 
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correlation distances along each of the orthogonal axes are significantly 
different. 

The correlator output where the template was obtained by convolv- 
ing the grid operator G with is shown in Fig. 4-12(c). The peak 
of the correlation is readily identified as in the previous figure; 
however, this processor exhibits a significant suppression of the 
pedestal variance. 

An urban scene is shown in the outlined area of Fig. 4-13(a) 
and the corresponding correlator outputs are shown in Figs. 4-13(b)-(e). 
As was previously observed for the agricultural area, the location 
of the peak in the correlator output using the image template is 
readily established. However, the pedestal variance in the correlator 
output using the processor template is significantly reduced. 

The outlined area in Fig. 4-14(a) is an example of a natural area. 
The correlator outputs shown in Figs. 4-14(b) and (d) are significant 
because in each of the outputs the peak corresponding to the correct 
juxtaposition of S on the background appears as a narrow peak project- 
ing from a larger bump. If a hill-climbing algorithm was being used 
to locate the correct juxtaposition of S and the step size was large, 
it is conceivable that an incorrect maximum would be found. However, 
in Figs. 4-l4(d) and (e) these undesirable lobes have been reduced 
and the central peak is readily identified. 

The significant numerical quantities associated with these results 


are tabulated in Table 4-4. 



Figure 4-13 Correlator Output for an Urban Area 




(d) Image Template 

N 2 - 4096 


(e) Processor Template 

N 2 = 4096 



Figure 4-13 Cont . 





(b) Image Template 
N 2 - 1024 


Figure 4-14 Correlator Output for a Natural Area 




(d ) Image Template 
N 2 = 4096 


(e) Processor Template 
N 2 - 4096 



Figure 4-14 Cont. 
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Table 4-4. Significant Numerical Quantities in the 
Noise Free Line Scanner Data. 

Figure Template Adjacent Cell Signal- to- 

Size Correlation Noise Ratio 


4- 12(a) 

64 x 64 

4-12(c) 

64 x 64 

4-13(b) 

32 x 32 

4-13(c) 

32 x 32 

4-13(d) 

64 x 64 

4-13(e) 

64 x 64 

4-14(b) 

32 x 32 

4- 14(c) 

32 x 32 

4- 14(d) 

64 x 64 

4-14(e) 

64 x 64 


.775 

133 

.775 

4510 

.629 

106 

.629 

664 

.690 

234 

.690 

1880 

.577 

21.8 

.577 

75.9 

.570 

10.4 

.570 

966 
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4.3 Multitemporal Data 

Multitemporal data sets generally exhibit both misregistration 
and noise components. The processor template, given by equation (3-78), 
requires both the signal and noise covariance matrices. In the general 
case, these matrices must be estimated from the data, although for 
certain classes of data the form of the matrix is assumed and the 
parameters needed for determining the matrix elements are estimated 
from the data. 

The difficulties in obtaining estimates of these covariance 
matrices were discussed in Chapter 3. Briefly, however, these diffi- 
culties arise from the requirements of large arrays of computer memory 
and the size of the data sets necessary to obtain meaningful statistical 
estimates of these matrices. 

In view of these computational difficulties the development of the 
template given by equation (3-78) will not be carried out; instead, the 
processor used previously in this chapter will be used to experimentally 
evaluate its applicability to multitemporal data sets. 

The data sets selected for multitemporal processing are shown in 
Fig. 4-15 with the reference data S selected from set A and the back- 
ground data £5^ chosen from set B. Four areas are selected from this 
data for processing, and the approximate locations of these areas are 
outlined in data set A. The following procedure is used in analyzing 
these data sets. The correlator outputs for these test areas for 
values of the prewhitening coefficient of 0.0 and 0.9 are computed. 

From these results it is concluded that the algorithm will not give 
meaningful results for geometrically uncorrected data regions. 
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Table 4-5. Characteristics of the Experimental 
Multi temporal Data 


Purdue Flight Line 210 



DATA SET A 

DATA SET B 

Run No. 

71053900 

71062900 

Date 

13 Aug 1971 

30 Aug 1971 

Time 

1202 

1229 hours 

Altitude 

5000 

5000 feet 

Heading 

180 

180 deg. 

No. Samples/Line 

222 

222 

Angular Resolution 

6.06xl0" 3 

5.82x10 3 radian 

Wavelength Band Used 

.54 - .60 

.54 - .60 ^ 
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The moduli of the Fourier transforms of these data regions are 
then computed, and the data of Region 2 geometrically corrected. The 
section concludes with the crosscorrelator outputs where the corrected 
data are the crosscorrelator inputs. 

4.3.1 Crosscorrelation of Misregistered Data 

The outputs of the crosscorrelator for the data sets selected 
from the set of multitemporal imagery are shown in Fig. 4-16. For the 
results shown in this figure, the leftmost display is for a value of 
the prewhitening coefficient a - 0.0. This corresponds to using the 
reference data set as the correlator template. The outputs for 
cv =» 0.9, a value chosen to approximate a nominally correct value of 
this coefficient for most previously observed data, are shown in the 
rightmost pictures. 

The correlator outputs shown in Fig. 4-16 for regions 2 and 3 
exhibit the typtically observed results for scenes consisting of 
rectangular fields. Observing the location of region 1 in Fig. 4-15, 
it is seen that this area too is composed of agricultural fields. 
However, data set B exhibits a non-linear distortion component in 
this area, suggesting that the aircraft began executing a maneuver 
during the period of time in which these data were being taken. 

The correlator output shown in Fig. 4-16, region 4, shows a 
results obtained with input data selected from an area containing 
undeveloped terrain. This output exhibits a large, readily distinguish- 
able unimodal peak. With data from a class exhibiting such a cross- 
correlation function, efficient algorithms can be developed such that 




Figure 4-1 6 Correlator Outputs for Multitemporal Data 








Figure 4-1 6 Cont . 
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a complete search of a region of the correlation output plane for 
a maximum is not required. Rather a search procedure is implemented 
where successive steps move toward the peak. 

The results for a - 0.9 without exception are of little value. 
However, for the processor to yield useful results, it will be demon- 
strated that it is necessary to implement a geometrical correction of 
the data. 

4.3.2 Regional Geometric Distortion 

Regional geometrical distortion can, for certain classes of dis- 
tortion and data, be estimated by determining the differences in the 
moduli of the two-dimensional Fourier transforms of the two approximately 
corresponding regions under consideration. Examples of moduli of 
Fourier transforms of line scanner data are shown in Fig. 4-17. The 
Fast Fourier Transform algorithm was used for computing the transform 
of these data. In each of these pictures the transformed data have a 
logarithmic amplitude scale, and the gray scale employed in the display 
system is linear. The size of the data set being transformed is 
128 x 128 picture elements. 

The transforms of regions 1 and 3 exhibit the concentration of 
spatial- frequency energy in the linear loci typical of structured data 
scenes, whereas the moduli of the transforms of region 4 exhibit 
unstructured characteristics. In each case it is of significance that 
most of the energy appears to be concentrated in the very low 
frequencies. The double lines in the moduli of the structured data 
are due to the aliased energy from the pattern located at the adjacent 
harmonics of the sampling frequency. 


Data Set A 
Region 1 


Data Set B 
Region 1 


Data Set A 
Region 3 


Data Set B 
Region 3 


Data Set A 
Region 4 


Data Set B 
Region 4 


Figure 4-17 Moduli of Two-Dimensional Fourier Transforms 
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The data from region 2 were selected as the set to be geometri- 
cally corrected. The moduli of the transforms of these data are 
shown in Fig. 4-18. The algorithm for locating the loci of energy 
in the transform domain, which was developed in Chapter 2, was 
applied to these transformed data, and the angular quantities given 
by the algorithm, measured with respect to the u^ axis, are as follows 

-2 

data set A -7.18 x 10 radian 

1.57 

data set B -3.30 x 10 ^ " 

1.52 


With a hypothetical transform coincident with the coordinate axes 
used as a reference data transform, the resultant angular quantities 
defined in Fig. 2-13 are 


data set A 


data set B 


» -4.1° 

0 2 = 0 ° 

- -1.9° 
0 2 - -3.3° 


The two-dimensional polynomial given in equation (2-14) was 
used to implement the correction. Rather than mapping one data set in- 
to the other, it was chosen to rectify each data set assuming the 
reference data set was coincident with the coordinate axes. Examining 
the data sets in both the spatial and spatial- frequency domain, the 
angular quantities chosen for the rectification differed slightly 
from those above. The numerical values used for the rectification are 



Data Set A Data Set B 


Figure U-l8 Modulus of Distorted Data and Geometrically 
Corrected Data for Region 2 
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data set A 


data set B 


with the reference transform angular quantities having the values 
9 1 - 0 ° 

0 2 = 90° 


for each of the data sets. 

Thus for data set A using equation (2-56) the value of the skew 
coefficient is 


t,n [ 2 180 


- to "[lSo (3 ' 7 >]- 


The distortion matrix, defined by equation (2-58), is then 


0.9979 -0.0645 


1.002 


The corresponding quantities for data set B are 


tan [iio <86 - 7 + 1 ' 7) ] 


tan [ 180 (1 - 6) ] “ °- 0279 


and the distortion matrix is 
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-0.0297 
0.9987 

The results of this geometrical correction on the spatial data 
are shown in Fig. 4-18. It is emphasized that the correction is 
applicable only in the center of each of these pictures as the 
distortion introduced by the aircraft motion is a dynamic quantity. 

4.3.3 Correlator Output for Geometrically Corrected Data 

The output of the crosscorrelator for the corrected data set is 
illustrated in Fig. 4-19 for various values of the prewhitening coef- 
ficient O'. The template size for each of these outputs is 64 x 64 
picture elements. Comparing the result of Fig. 4-19(a) to the corres- 
ponding output of region 2 of Fig. 4-16, it is seen that the geometrical 
correction has increased the sharpness of the peak at the correction 
justaposition of 

A significant result is that the pedestal variance is reduced while 
the central peak becomes more discriminable as the value of a approaches 
unity. A plot of the correlator output signal-to-noise ratio for these 
results is shown in Fig. 4-20. 


—2 


0.9999 


0.0576 



Figure 4-19 Correlator Output for Corrected Multitemporal Data 






Figure 1+-19 Cont . 
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Figure 1^-20 Correlator Outpu 
Function of the 
(Geometrically C< 
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CHAPTER 5 
SUMMARY 

5.0 Summary 

Given two sets of imagery of some common area taken at different 
times, it is observed that these data exhibit geometrical misregistra- 
tion. The rectification of this imagery requires that the parameters 
be determined for the mathematical model assumed to characterize the 
misregistration. The problem studied in this investigation is the 
identification of conjugate, or corresponding, points in multi- 
temporal imagery such that the coefficients of the misregistration 
model can be obtained. 

Selecting a reference data subimage, denoted S^, from the first 

data set, these conjugate points are located by finding the correspond 

ing subimage S/ which is embedded in the background data S^. The 

measure of similarity between data sets is taken to be amplitude 

correlation. Assuming the misregistration between S and S* to be 

— r — r 

modeled by an affine transformation, a sequential procedure is 
implemented in which the geometrical distortion components are first 
estimated and removed. Algorithms utilizing the moduli of the two- 
dimensional Fourier transforms of the data S^, and its approximate 
corresponding subimage are developed for estimating the distortion 
model coefficients. 



Following the geometrical correction, the displacement 
coefficients are estimated by determining the location of the peak 
of the correlation surface computed by crosscorrelating a template, 
generated from the reference data, with the background S^. The peak 
of this correlation surface is assumed to be the point of correct 
superposition of the reference data on the background data. 

The algorithms reported in the literature for computing the 
correlation surface, or some equivalent quantity, do not give a 
surface on which the point of greatest correlation is maximally 
discriminable from the values at other points of the surface where 
the criterion of discrimination between points on the surface is as 
defined in Chapter 3. The one exception is the filter reported by 
Arcese jet al. [6] . 

The filter derived in this investigation extends these previous 
results to data having a covariance matrix K g and in which the back- 
ground data contains an additive noise component with statistical 
properties characterized by the covariance matrix K^. The analytical 
expression of this filter is obtained under the constraints of both 
an output of unity at the correct superposition of the reference data 
on the background and for a bounded energy in the filter. With the 
removal of both constraints and for no additive noise components, this 
filter reduces to the result of Arcese ejt til. If in addition the 
correlation between adjacent picture elements is neglected, which 
yields the identity matrix for K g , the filter then becomes the 
reference image itself. 
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The application of this filter to noise free and distortionless 

data, with the covariance matrix of the image assumed to be of a 

simple exponential form, yields results which substantiate the 

theoretical development of the processor. In these results the 

correlation surfaces of the data sets processed with the filter 

exhibit narrow central peaks and have significantly decreased variance 

in the pedestal region of the surface. 

The application of this algorithm to multi temporal data 

introduces some problems which at the present have no solution. 

Implementation of the general filter for this data requires knowledge 

of the covariance matrices of both the reference and noise processes; 

however, it is not possible to make meaningful estimates of these 

quantities at the present time. The difficulty in estimating these 

matrices arises from two sources. The first is the size of the arrays 

required for handling the matrices in the computer. It has been 

experimentally found that the useful sizes of reference data sets 

range between 2^ x 2^ to 2^ x 2^ picture elements, with the resulting 

sizes of the covariance matrices ranging between 2^ x 2^ and 
14 14 

2 x 2 elements. The second difficulty is the unavailability of 

a sufficiently large set of data having similar statistical properties 

to the reference such that a meaningful estimate of the covariance 

matrix of the reference data can be made. An additional difficulty 

in estimating the noise properties is that, with the noise defined 

as the difference between the data sets, the determination of K 

~n 

requires that the data be in registration to obtain the differences 
which are related to scene changes and are not due to misregistration 


error. 
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If It is possible to approximate the estimated covariance matrices 
with matrices having a Topi itz form, the requirement for computer 
storage can be reduced by a factor of l/m where the number of elements 
in the reference data set is m x m. However, the requirement remains 
for the availability of data sets of sufficient size from which to 
make the initial covariance estimates. 

The use of the algorithm derived for the noise free processor 
with both noise free and multitemporal data yields encouraging results 
in that the central peak of the correlator output Indeed becomes more 
discriminable from the output in the pedestal region where discrimina- 
tion is defined as the ratio of the correlator output signal-to-noise 
using the optimum filter to the corresponding correlator output 
quantity with the reference scene used as the filter. The increase in 
discrimination for the noise free data ranges between 3 to 15 dB and 
for the multitemporal data the corresponding increase is greater than 
6 dB. 

A requirement imposed on the data by the filter for the optimum 
discrimination of the correct superposition of is that the 
geometrical distortion components of this misregistration model must 
be reduced below the magnitudes found to typically occur in the multi- 
temporal data. It is demonstrated in this investigation that the 
linear geometrical distortion existing between two data sets, can be de- 
termined from the moduli of the Fourier transforms of these data. Values 
of the skew and rotational distortion components were determined for 
a set of multi temporal imagery using the linear regression algorithms 


155 


developed In Chapter 2. These values of these distortion components 
were found to agree with numerical values obtained by independent 
methods to within a few tenths of a degree. 
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APPENDIX A 

THE TWO-DIMENSIONAL FOURIER TRANSFORM 

The properties of the two-dimensional Fourier transform of the 
function f(x) under an affine transformation 

H-Ax+t (A-l) 

where A is a [2 x 2] nonsingular matrix and 



are of interest in this study. In the following, expressions for 
both cartesian and polar coordinate systems are derived and the 
symmetry between the spatial and spatial- frequency domains is 
discussed. 

The two-dimensional Fourier transform of the function f (x) is 
defined as 

CD 

F(u) ■ J J f (*> ex P ["J2 tt (u,x) 1 dx (A-2) 

-CO 

and the inverse is 

f (x) - J J F (“) exp T j2n (u,x) 1 du (A-3) 

—00 

where 



162 


is the spatial frequency and (u,x) is the Inner product 

T 

(u,x) ■ u x • Uj Xj + Uj * 2 ' . 

(a) The Affine Mapping 

The Fourier Transform of the function 
H 


f(x) 

* f (H X) - g(2) 

(A-4) 

F(u) 

CO 

- J J f(H x) exp ^-j2n(u,x)^J dx 

— CO 

(A- 5) 


“ "IjJ J J ex P {"J 2rr [(H.a" 1 ! - A* l t)^j 

} di (A-6) 


“ fjf exp [ a" 1 t/j 



• J J exp » a’ 1 i>] d i 

— CO 

(A-7) 


m Tjy exp [ -J 277 (h» a" 1 t)J 



* J I 8 <£> ex P [(A* 1 u, x)] } A X (A- 8) 


» jjy exp [ J2tt (u,t)^j ? [(A _1 t) uj (A-9) 

- G(v) 

where J is the Jacobian of the transformation 


tyl dy 



^y 2 ^y 2 

dx^ ■ b-x.2 


(A- 10) 
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The relation between these coordinate variables Is 
f(x) < > g(i) 

I ft 

F(u) < > G(v) 

where denotes a Fourier transform pair. 

(b) The Two-Dimensional Fourier Transform in Polar Coordinates 
To obtain the equation analogous to equation (A- 10) in polar 
coordinates, let 

Xj » r cos 9 (A- 11) 

Xj ■ r sin 9 

and 

u^ « p cos 0 

(A- 12) 

u 2 - p sin 0 

Substituting equation (A- 12) Into equation (A-2) 

«» 2tt 

F(u) - J J r f[x(r,9)]exp £-j2n (u,x)J d0 dr (A- 13) 

o o 

» 2rr 

"II r ex P £”j2n( u i r cos ® + u 2 r 8in ®)J d ® dr 

o o 

(A- 14) 


Now using equation (A-ll) in (A-14) 


F[u(p,0)] »G(p,0) 

00 2tt 

- I I r g(r,9) exp ^-j2nrp cos (9-0) d9 dr (A-15) 

o o 
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(c) Scale and Rotation In Polar Coordinates 

Uniform scale changes or rotations of the function g(r,9) are 
related by the coordinate mapping. 

T 

g(r,0) -♦ g(ar,6 + e Q ) (A-16) 

Substituting equation (A-16) into equation (A- 15) yields 

OD 2tr 

G'(p,0) ” J J r g(ar, 9 + 9 q ) exp £-j2rrrp cos(0-0)J dO dr (A-17) 

o o 

Letting £ “ ar 


X - 9 + 0 

__ o 


» 2rr ° 

G # (p,0) - J J (1) g(S,X) exp | 

-J2T, (S)p 

cos (X-e o -0)]dX d0 (A- 18) 

o o 



“ 2tt 



” 7“2 J ] 5 8 (5 > x) 

a 2 J J 
1 1 o o 

p(-J2n5 (£ 

| cos [x - C0 +0 o )]j} dXd5 


(A-19) 

- G(p/a, 0 + 8 > 

1*1 


(A- 20) 


(d) Spatial and Spatial-Frequency Domain Symmetries 

For the affine transformation defined by equation (A-l) the 
resultant transform was given by equation (A-9) as 

G(r) ■' Jjj; exp £ J2rr(u,t)] ? [(A _1 t) u| (A- 21) 

Thus the coordinates are related as 

v » (A -1 ) 1 u (A- 22) 
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The energy density spectrum of f(x) is 

S(u) - |F(u)| 2 (A-23) 

and the linear phase term cancels. Thus the coordinates of the energy 
density spectrum (or equivalently the modulus) are related to the 
linear distortion matrix in the spatial domain, whereas any spatial 
shift results in a linear phase term in the spatial frequency domain. 
For real signals it follows immediately from equation (A-2) that 

F(-u) - F*(u) (A-24) 

where the asterisk denotes the complex conjugate. Thus the energy 
density spectrum (or modulus) is symmetric about any line passing 
through the origin in the spatial frequency domain. 
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APPENDIX B 

FILTER NORMALIZATION COEFFICIENTS 

The normalization coefficients for the filters used to smooth 
the random number ideal data are derived assuming a continuous 
functional form. These weighting coefficients are then used with 
the discrete two-dimensional convolution expression to implement 
the smoothing filter. As a check for the correctness of these 
coefficients, the gain (volume) of each filter was experimentally 
verified by convolving each filter impulse response with an array, 
each element of which was unity. 

Separable Exponential Filter 

The filter impulse response is given by 

h^x) - k x exp ~ ^ + *2>J ( B_1 ) 

R 0 > °» |*il. I x 2 i < R < " ♦ 

The gain of h^(x) is 
R R 

G i ■ 4 JI k i eicp [- 1; (x i + * 2 >] d *i dx 2 <B ' 2 * ) 

o o 

■ 4 k i R o[ l * ** (^) ] 


(B-2b) 
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For unity gain 

k i * { 2 R o[* - « xp (r)]} 

Isotropic Exponential Filter 

The isotropic filter impulse response is given by 

h 2 (x) - k 2 exp (- ) 

o 

2 2 1/2 

r - <x^ + x 2 ) > 0, 0 < R q < R < °°. 

The gain of h 2 (x) is then given by 
R 2tt 

®2- J J k 2 r “P (- r ) d9 dr 

o 

O o 

R 

- 2tt k 2 J r exp (- f- ) dr 

o 

o 

let 


X - f- , R dX - dr 
o 

Thus 

R / R o 

G 2 » 2tt k 2 J R^ X exp (-X) dX 


and integrating by parts 

R/R ®"/®o 

G 2 -*» 2tt k 2 R^ [-X exp (-X) \ ° + J exp (-X) dxj 

o 

o 

- * k 2 R o [- r ex " (- f ) - (- f ) + 

o o o 

- 2f rR k 2 [l - (l ♦ f ) .xp (- f ) j 


(B-3) 

(B-4) 

(B-5a) 

(B-5b) 

(B-6) 

(B-7) 

(B-8a) 

(B-8b) 

(B-8c) 


■ 168. 

For unity gain 

k 2- {* R o b - b * r ) {■¥)}' »- 9 > 

o o 

Test of Gain 

For each filter the convolution 

C(u, v)' « E I h(u-i, v-j) a(i,j) (B-1’0) 

i j, i u,v*»0 

was computed with the filter having characteristic length R q « 1.0 and 
an a(i,j) consisting of a 100 x 100 array of element each element, 
having a value of unity. The values obtained are 

kj - 1.0006 

k 2 - 0.9998 . 
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APPENDIX C 

CORRELATOR OUTPUT SIGNAL TO NOISE RATIO 


In the absence of uncorrelated additive noise the correlator 
output signal- to-noise ratio was shown in Chapter 3 to be given by 


| - S T K~ 1 S 
N — r — — r 


(C-l) 


where S is the reference signal and K is the covariance matrix of 
the unit variance process of which ^ is a finite realization. It 
is assumed that the covariance matrix has the form 

.n-1. 1 


where 


and 


K - 


k 

pk 


pk 

k 


p 2 k 


Pk 


n-l. 

p k • • • 


k- II p tJ I I 


f 1 , i-J 

'^^Nl , i^. 


(C-2) 


(C-3) 


(C-4) 


where p is the correlation between adjacent cells. 
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Hie covariance matrix then has the explicit Inverse 


K 


•l jl_ 

i-p : 


-l 


•P.k 

0 


-1 


-Pi - 1 

(i + p 2 >k' 1 


-pk 


-1 


-Pk 

(1 + p 2 )^’ 1 


(C-5) 


where 


, -1 1 
k - ; 

lV 


1 

-P 

0 


-P 0 

2 


1 + p 
-P 


-P 

1 + P J 


1 J 


(C-6) 


To evaluate equation (C-l) , substitute equation (C-5) into (C-l) 
and expand, noting that the r subscript has been deleted. 


S 1 
N " i 3 

1 -p 


(S T (l),...,S T (m)) 

If 1 

-Pk " 1 

0 

0 

S( 1 ) 


-Pk " 1 

(l+p 2 )k _1 

-Pk " 1 


• 

e 


o 

•• 

-Pk ” 1 

( 14 p 2 )k “ 1 

• 


S(») 


• 

• 

_ 0 



• 

S’* 



(C-7) 

-^~2 (S T (l)]f 1 - pS T ( 2 )k‘ 1 , -p/cDk ” 1 + (l+p 2 )S T ( 2 )k - 1 - p,S T ( 3 )k” 1 , 
1-P 


■pS^m-Dk" 1 + S T (m)k _1 ) 


S(l) 

S(m) 


(C-8) 


• • • » 
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| - ~ o Fs T (l)i “ 1 S(l) - pS T (2)k“ 1 S(l) - pS^Djf 1 S(2) 

1-p L 

+ (1+p 2 ) -S T (2)*f 1 S(2) - pS T (3)k _1 S(2) - pS T (2)k' 1 S(3) 

+ (1+p 2 ) < S T (3)Jt“ 1 S(3) - pS T (4)k _1 S(3) - 

... - pS^m-Dk” 1 S(m) + S(m)] (C-9) 


Now each term is a scalar, thus 


S T (i)k _1 S(J) - [S A (l)k _A S(j)] A - S i (J)k*’ A S(i) (C-10) 

since (k ■ k k ^ symmetric. 

| [SW 5(1) ' 2 P£^ (l)k~^ S(2) + p Z S X (2)k" 1 S(2) 

1-P 

+ S T (2)k _1 S(2) + ... + S T (m)k" 1 S(m)] (C-ll) 


T/ j « / tV 


f I S(i)k -1 S(i) - 2p I S(i)k _1 S(i+1) 
1-p L i-1 i-1 

m-2 T T "1 

I S 1 (i)k 1 ' S(i) 
i-1 J 


+ P 


Expanding a term of equation (C-12) 

S T (i)k -1 S(i) - -i-5 (s ,,...,8 n ) 
1-p 


(C-12) 


“P 

0 


-P 

l 

-P 


0 

-p 

1 


’il 


in 


1 


as !w 
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“ 2 [®il ' pS l2 V ‘ P8 11 


+ (1 + p ) b 12 - pa l3 , 


2 X • 


’ P8 i2 + (1 + P ‘> 8 i3 ' p8 i4’*"*’‘ p8 i,n-l + 8 in] 


~ L 2 \*ll " 2 P: , ^i“i9 + (1 + P 2 )8 ?o - 2p s 4 ,s 

1-p L 


'il"i2 


+ (1 + P ) 


2. 2 


’13 


’i2 "12 13 

, 2 

- p 8. ,8. +8. 

v 1 ,n-l in In 


_ 


11 


in 


(C-14) 


(C-15) 


1 r N 2 N_1 2 N “ 2 2 I 

-SI l8 iJ ' 2P 1 8 tJ 8 i j+l + p S 8 iJ J 

1-p L j»i 1 *J TA j-1 J 


(C-16) 


[n - 2p(N-l) + p 2 (N-2)] 


1-P 


M 


(C-17) 


(C-18) 


N 


where N (— » ... 

j-1 J 


2 s 2 ) - N, 


since 5> is assumed to be from a unit variance 
N-l 


process and (N-l) (^[ ^ 8 ij 8 i j+l) “ P CN-1) where p is 


the 


correlation between adjacent cells. 

Substituting equation (C-18) into equation (C-12) yields 


M 


s i r “ M_1 2 M - 2 l 

£ - I N - 2p E pN + p* I N 

1-p L i-1 i-1 i-1 J 


(C-19) 
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Thus the processing gain of the two-dimensional crosscorrelator 
theoretically is equal to the total number of samples processed and 
is independent of the adjacent cell correlation p. 



APPENDIX D 

THE MAXIMUM PRODUCT OF A CONSTRAINED SUM 


If the sum of a set of real numbers X^ is constrained to 


the fixed value 


£ X ** K 


m 


m»l 

then the product 


P - H X. 
n*»l 


m 


is maximized for all X^ equal. 


Solving for X 


j 


n 


X. ■■ K - £ X_ 

J B-1 " 


n 


(k - £ X ) X. 

V wl *' 1 


* X j-1 X j+1 •** X n 


K(X ^ ... X j_j^ X j+1 * * * ^ n ) 


dP 


n 


■ J X (X . • • # X * ■» X « * « ••• X ) 
n.1 m 1 J' 1 J +1 n 




n 


ax 


j+i 


K(X 1 ... X m (X l*** X j-l X j+2 , *' X n ) 

“ ^ X l** ,X j-l X j+l*** X n^ 
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have 

(D-l) 

(D-2) 

(D-3) 

(D-4) 

(D-5) 

(D-6) 

a* 0 
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But 


K 


E X_ 

ra*l 


m 


therefore , 


EX - £ X - X... 
m-1 m „■! ° ^ +1 


- 0 


J^o 


(D~7) 


X j ’ X J+1 “ 0 


X J ■" X j+1 


(D-8) 

(D-9) 


for all j ■ 1, n-1. 

Thus all X t must be equal and have the value 


K 

n 


(D-10) 


APPENDIX E 
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THE AFFINE TRANSFORMATION 


The affine transformation Is defined as 
j^-Ax + jt-Bx 

where for a real two-dimensional space 



and 


a il a l2 
. a 21 ®22 


(E-l) 


Is a non-singular matrix. 

Two properties of this transformation are of importance in this 
study; these are proved in the following. 

(1) Any non-singular affine transformation maps lines into lines. 

Proof: In any affine space the vector x lying on the line 

passing through u, v can be represented as 

x ~ (1 - X) u + X v, X real (E-2) 

Then 
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B x ■ A [^(1 - X) u + X v J + t 

» (1 - X) (A u + t) + X (A v + t) 


(E-3) 


*»(l^X)Bu + XBv 

Thus B maps the line passing through u,v into the line passing 
through Bu, By . 

(2) Ratios of lengths of line segments are invariant under an affine 
transformtlonl 

Proof: u. u , v. v are elements in the real n-dimensional 

Euclidean space E n * The distance between vectors is defined as 



Consider the mapping 
B 

u * B u (E-5) 


with 

Bu-Au+t 

where A is an n x n nonsingular matrix and J: is the n x 1 displacement 
vector. 

Now the ratio of line segments in the affine space is 



R - 


!lH(u - u )|| 




(E-6) 


1|A u + t - A u - t|| llA(u - a >11 



where 1|A|| is the norm of the operator A. 
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